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ABSTRACT 


A  pool  of  liquid  with  a  horizontal  free  surface  is  bounded  on  one  side  by  a 
vertical  solid  wall,  which  is  maintained  at  a  cold  temperature  relative  to  the  core  flow 
region.  Strong  temperature  gradients  along  the  surface  give  rise  to  surface  tension 
variations  (thermocapillary  stress),  which  drives  flow.  Thin  viscous  boundary  layers 
form  along  the  surface  and  wall.  A  boundary-layer  model  is  designed  which  captures 
the  dynamics  of  the  cold  corner,  applicable  for  any  Marangoni  number  M  and  Prandtl 
number  P  in  the  convective  inertial  regime. 

Analytical  expressions  for  the  velocity  and  boundary-layer  thicknesses  are  de¬ 
veloped,  which  allow  accurate  prediction  of  the  flow  field.  The  core  flow  region  (out¬ 
side  the  viscous  boundary  layers)  is  treated  as  irrotational  flow  and  Laplace’s  equation 
is  solved  using  both  a  Green’s  function  approach  and  a  complex  variables  approach 
in  the  quarter-plane.  The  flow  along  the  wall  is  treated  as  a  plane  wall  jet. 

The  two-dimensional  unsteady  heat  equation  is  solved  using  an  alternating 
direction  implicit  method.  Results  show  that  the  flow  into  the  corner  is  strong  enough 
to  contain  the  thermal  field,  compressing  the  isotherms  along  the  wall  after  steady- 
state  is  reached.  Additionally,  a  uniform  stream  function  prediction  is  developed,  by 
matching  the  inner  and  outer  flows,  giving  a  relatively  accurate  depiction  of  the  flow. 
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I. 


INTRODUCTION 


One  must  learn  by  doing  the  thing;  for  though  you  think  you  know  it,  you 

HAVE  NO  CERTAINTY  UNTIL  YOU  TRY.  —  Sophocles 


The  purpose  of  this  research  is  to  examine  the  velocity  and  temperature  fields 
in  a  two-dimensional  fluid  flow  problem  embodying  a  thermocapillary  feedback  mech¬ 
anism.  The  structure  of  the  velocity  and  thermal  fields  will  change,  based  upon 
whether  the  heat  transfer  is  conductive  or  convective  in  nature,  as  well  as  whether 
inertia  plays  a  dominant  role  or  not.  One  of  the  governing  parameters,  the  Marangoni 
number  M,  measures  the  importance  of  thermal  convection  relative  to  thermal  dif¬ 
fusion.  A  second  governing  parameter,  the  Prandtl  number  P,  gives  the  ratio  of 
viscous  to  thermal  diffusion  for  the  material.  In  this  research,  we  study  the  effects 
of  a  high  Marangoni  number,  indicative  of  convective  heat  transfer,  combined  with 
a  low  Prandtl  number,  indicative  of  inertial  flow.  The  most  important  regime  for 
materials  processing  is  that  regime  where  thermal  convection  is  dominant  and  the 
flow  is  inertial,  with  viscous  boundary  layers  forming.  This  dissertation  outlines  a 
model  for  this  important  regime,  with  the  goal  of  predicting  the  velocity  and  thermal 
fields  which  will  occur. 

Chapter  II  gives  a  background  on  the  importance  of  this  research  and  possible 
physical  applications,  followed  by  a  survey  of  related  work  in  thermocapillary  flows  by 
previous  researchers.  In  our  literature  search,  we  did  not  find  any  other  researchers 
who  had  studied  the  cold  comer  region  with  large  Marangoni  numbers  and  small 
Prandtl  numbers.  A  representative  velocity  vector  field  (using  Canright’s  numerical 
data  [Ref.  1]),  which  serves  as  a  guiding  principle  for  our  assumptions,  is  shown  at 
the  chapter’s  end. 

Chapter  III  formulates  the  problem  mathematically  and  gives  all  assumptions 
and  a  thorough  scaling  analysis.  It  is  important  to  obtain  a  summary  of  theoretical 
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scales  of  the  thermal  length,  viscous  thickness  and  length,  and  surface  velocity.  Al¬ 
though  the  scalings  derived  are  different  than  those  previously  published,  they  agree 
better  with  the  published  numerical  data. 

Chapter  IV  contains  the  conservation  equations  and  how  they  affect  each  of 
four  regions  within  the  cold  corner.  Using  the  appropriate  scaling  factors,  once  the 
velocity  is  determined  in  one  region,  it  is  matched  to  that  in  all  neighboring  regions. 

Chapter  V  discusses  the  approximation  methods  employed  to  describe  the 
velocity  profiles  in  the  viscous  boundary  layers  and  core  flow  region.  Two  methods 
from  the  literature  were  modified  to  predict  velocity  behavior  along  the  free  surface 
and  the  rigid  wall.  The  modification  along  the  free  surface  introduces  a  new  technique 
to  predict  the  velocity  from  a  free  surface  with  a  thermal  gradient  boundary  condition 
and  no  non-slip  condition  through  a  viscous  boundary  layer  into  the  core  flow  region. 
The  flow  from  the  corner  down  the  rigid  wall  is  treated  as  a  plane  wall  jet.  An 
expression  for  the  temperature  solution  in  the  domain  is  then  determined. 

Chapter  VI  shows  the  numerical  results  obtained  from  various  techniques. 
The  alternating  direction  implicit  (ADI)  method  is  the  basis  for  solving  the  two- 
dimensional  unsteady  heat  equation.  Various  tests  against  problems  with  known 
solutions  are  performed  to  validate  the  code,  followed  by  the  solution  of  the  present 
problem.  The  scheme  is  also  modified  to  update  the  velocity  profile  during  the  time- 
step  iterations  to  model  the  physical  feedback.  A  modification  to  a  wall  jet  profile  is 
incorporated  into  the  scheme.  The  results  are  presented  as  a  uniform  velocity  field, 
combining  the  boundary  layers  and  core  flow,  for  comparison  with  previous  numerical 
results. 

Chapter  VII  gives  the  conclusions  and  a  few  discussions  associated  with  them. 
Chapter  VIII  fists  selected  areas  of  possible  future  research. 

Appendices  A,  C,  and  D  show  work  which  attempted  to  improve  the  model  in 
the  surface  and  wall  boundary  layers.  We  ultimately  did  not  use  these  methods  (for 
reasons  given  in  each  appendix),  but  we  felt  we  should  include  them.  Appendix  B 
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gives  a  derivation  of  an  important  and  often-used  constant.  The  numerical  codes  are 
listed  in  Appendix  E. 
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II.  BACKGROUND 


We  do  not  set  out  on  this  path  of  enquiry  just  as  adventurers.  We  follow 

THE  FOOTSTEPS  OF  SUCH  GIANTS  AS  EULER  AND  LAGRANGE,  WHO  ALSO  WONDERED 
WHERE  IT  LED,  AND  PRANDTL,  WHO  WALKED  IT  AND  FOUND  THE  ANSWER.  —  David 
Pnueli  and  Chaim  Gut  finger 


A.  MOTIVATION 

Materials  processing  often  involves  the  melting  and  solidification  of  the  mater¬ 
ial.  Several  practical  processes,  such  as  welding,  containerless  processing  of  materials, 
float-zone  purification,  and  Czochralski  crystal  growth,  involve  a  pool  of  molten  metal 
with  a  free  surface  [Ref.  1].  The  purest  large- volume  silicon  single  crystals  are  being 
commercially  grown  by  floating-zone  melting  [Ref.  2,  3].  Inevitably,  there  exists  a 
fluid  flow  within  this  molten  region. 

When  heavier  fluid  is  on  top  of  fighter  fluid,  the  weight  of  the  heavier,  denser 
fluid  makes  it  tend  to  sink,  while  the  fighter,  less  dense  fluid  tends  to  rise.  This  sets  up 
motion  within  the  fluid  known  as  buoyancy-driven  convection  or  natural  convection , 
also  known  as  Rayleigh  convection.  One  of  the  most  commercially  important  materials 
processes  that  are  subject  to  buoyancy-driven  convection  is  that  of  crystal  growth. 
A  major  advantage  of  performing  experiments  in  a  space  laboratory  is  being  able  to 
minimize  gravitational  effects  and  therefore  reduce  buoyancy-driven  convection  [Ref. 
4].  There  are  important  reasons  to  avoid  having  buoyancy-driven  convection  during 
crystal  growth.  Convection  can  alter  trace  element  distribution  in  the  resulting  crystal 
and  may  degrade  its  quality.  Many  other  physical  processes  can  greatly  affect  crystal 
growth,  such  as  radiation  heat  transfer  and  surface  tension  driven  flow.  Typically,  all 
of  these  processes  happen  at  the  same  time  on  Earth.  Buoyancy-driven  convection  is 
the  predominant  fluid  flow  on  Earth,  since  the  gravitational  forces  are  stronger  than 
surface  tension  forces  [Ref.  5].  In  microgravity,  however,  buoyancy-driven  convection 
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is  reduced  by  a  factor  of  10-6,  allowing  thermocapillary  flow  to  predominate  [Ref. 
6].  The  near  elimination  of  buoyancy-driven  convection  has  promise  for  materials 
processing  in  space,  a  low-gravity  environment.  For  such  applications,  it  is  essential 
that  thermocapillary  flow  be  studied  in  detail,  without  the  influence  of  buoyancy- 
driven  convection. 

Convection  in  the  molten  metal  is  typically  vigorous  and  can  greatly  affect 
the  results  of  the  materials  process,  including  the  size  and  shape  of  the  melt  pool, 
the  heat  transfer,  the  mixing  of  solutes,  and  the  final  microstructure  of  the  product 
[Ref.  7].  The  buoyancy  force  can  also  play  a  role  in  establishing  the  fluid  flow. 
When  the  dimension  of  the  molten  region  is  small,  as  in  the  case  of  welding  and 
laser  materials  processing,  the  buoyancy  force  becomes  unimportant.  Beyond  the 
molten  region,  the  absorbed  heat  is  conducted  into  the  substrate  and  transported 
away  by  the  motion  of  the  substrate.  Surface  tension  gradient  driven  flow  has  been 
identified  to  be  responsible  for  ripple  formation  in  the  weld  [Ref.  7].  When  buoyancy, 
electromagnetic,  and  surface  tension  forces  were  considered  in  a  numerical  solution 
for  the  convective  heat  transfer  within  the  molten  pool  during  arc  welding,  it  was 
found  that  the  surface  tension  gradient  is  the  dominant  factor  in  many  cases  [Ref. 
8].  Other  authors  have  verified  the  dominance  of  thermocapillary  flows  in  materials 
processing  [Ref.  9,  10,  11]. 

In  addition  to  the  heat  transport  due  to  the  motion  of  the  substrate,  there 
is  also  the  convection  due  to  the  fluid  flow  within  the  molten  pool.  Depending  on 
the  type  of  heat  source,  the  fluid  flow  may  be  caused  by  the  surface  tension  gradient 
[Ref.  7].  When  a  free  surface  is  heated  by  a  concentrated  heat  source,  the  resulting 
temperature  distribution  causes  a  nonuniform  surface  tension  distribution  [Ref.  12], 
Typically,  surface  tension  is  a  decreasing  function  of  temperature.  Thus,  the  fluid 
layer  on  the  surface  experiences  a  shear  force  pulling  the  fluid  from  the  warmer  central 
area  of  the  molten  pool  to  the  cooler  outer  region.  Since  the  bulk  fluids  are  viscous, 
they  are  dragged  along;  bulk  fluid  motion  then  results  from  the  temperature  gradients 
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along  the  interface.  This  phenomenon  is  known  as  thermocapillary  or  Marangoni 
convection.  The  Marangoni  effect  is  associated  with  movement  in  a  fluid  interface, 
caused  by  local  variations  in  interfacial  tension  that  are  caused  in  turn  by  differences 
in  composition  or  temperature.  A  drying  coat  of  paint  is  a  simple  example  where  local 
variations  of  surface  tension  set  into  intricate  motion  an  entire  liquid  film  [Ref.  13].  In 
the  case  of  the  float-zone  crystal  growth  process,  a  melt  zone  is  moved  slowly  through 
a  cylindrical,  vertically-oriented  crystal.  The  liquid  zone  between  the  polycrystalline 
incoming  material  and  the  single  crystal  that  is  grown  is  stabilized  by  surface  tension 
[Ref.  3] .  The  temperature  distribution  varies  along  the  interface  of  the  liquid  and  gas, 
which  leads  to  surface-tension  gradients,  creating  a  considerable  thermal  Marangoni 
flow  driving  force,  or  thermocapillary  convection. 

B.  PREVIOUS  RESEARCH 

There  have  been  many  theoretical  studies  of  thermocapillary  flows,  both  nu¬ 
merical  and  analytical.  Ostrach  [Ref.  5]  reviewed  of  some  of  the  analytical  studies  rel¬ 
evant  to  low-gravity  conditions.  Some  fluid  flows,  such  as  buoyancy-driven  and  g-jitter 
convection,  are  acceleration  dependent,  while  others,  such  as  surface-gradient,  ther¬ 
moacoustic,  and  phase-change  convection  are  independent  of  the  acceleration  field. 
There  are  forces  in  a  space  laboratory  which  act  on  experiments  in  a  gravity-like  man¬ 
ner.  G-jitter  (or  residual  acceleration  [Ref.  4])  is  the  term  used  when  combining  these 
gravity-like  forces  together.  Thermoacoustic  convection  refers  to  the  pressure-driven 
convection  that  results  when  a  confined  compressible  fluid  is  heated  rapidly;  the  name 
describes  the  sonic  character  of  induced  pressure  waves.  When  a  fluid  solidifies  or 
vaporizes,  a  change  in  density  associated  with  the  phase  transition  may  also  cause  a 
change  in  volume.  A  shrinkage  occurring  during  solidification  may  result  in  a  volume 
reduction,  further  resulting  in  the  flow  of  a  liquid  toward  the  solidifying  surface.  This 
is  known  as  phase-change  convection.  Davis  [Ref.  9]  reviewed  thermocapillary  in¬ 
stabilities  in  planar  layers.  He  gave  details  about  studies  concerning  Marangoni  and 
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hydrothermal  instabilities  involving  the  one-layer  system  geometry.  Kuhlmann  [Ref. 
3]  gave  a  thorough  reference  list  of  many  processes  subject  to  Marangoni  effects,  such 
as  the  motion  of  droplets  and  bubbles,  the  stabilization  of  thin  films,  evaporation 
and  boiling,  welding,  electrokinetic  effects,  and  crystal  growth  techniques,  to  name 
just  a  few.  Hondros  [Ref.  14]  offered  several  examples  of  phenomena  attributed  to 
Marangoni  convection  relevant  to  modern  technology,  such  as  in  hot  salt  corrosion  in 
turbine  blades,  the  drying  of  silicon  wafers  in  the  electronics  industry,  and  micro-pools 
produced  by  plasma  disruptions  in  a  prospective  thermonuclear  fusion  reactor. 

Cowley  and  Davis  [Ref.  15}  analyzed  the  two-dimensional  thermocapillary  flow 
near  a  hot  wall  for  vigorous,  viscous  flow  (large  Marangoni  and  Prandtl  numbers). 
The  fluid  flows  up  a  wall,  which  is  kept  at  a  constant  high  temperature  within  a  finite 
distance  from  the  free  surface,  and  then  turns  and  flows  along  the  free  surface.  This 
is  called  the  hot  corner  problem.  The  finite  distance  determines  the  length  scale.  The 
parameter  regimes  studied  either  have  no  viscous  boundary  layers  or  have  viscous 
boundary  layers  that  are  much  thicker  than  the  thermal  layers.  The  solution  they 
find  is  the  thermocapillary  analog  of  buoyancy-driven  convection  in  a  quarter  plane 
described  by  Roberts  [Ref.  16]. 

Sen  and  Davis  [Ref.  11]  studied  steady  flow  in  two-dimensional  slots,  which  are 
differentially  heated,  inducing  the  surface-tension  gradient  along  the  interface.  They 
studied  the  case  with  the  aspect  ratio  A  (ratio  of  the  depth  to  the  width)  approaching 
zero.  The  ends  of  the  slots  are  maintained  at  a  fixed  temperature  difference.  Flows 
on  the  interfaces  are  directed  from  the  hot  towards  the  cold  end  and  return  along 
a  region  removed  from  the  interfaces.  The  pressure  is  higher  at  the  cold  end  and 
the  interface  thus  bulges  near  the  cold  end  and  is  constricted  near  the  hot  end.  The 
leading-order  outer  solutions  having  parallel  flow  and  flat  interfaces  continue  to  be  the 
leading  order  outer  approximations  even  when  conditions  on  the  Reynolds  number 
and  Marangoni  number  are  relaxed  from  0(A). 

Zebib,  Homsy,  and  Meiburg  [Ref.  17]  conducted  numerical  studies  of  flow  in 
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a  two-dimensional  square  cavity  with  one  hot  wall  and  one  cold  wall.  They  showed 
that  for  moderate  to  small  Prandtl  numbers,  the  cold  corner  has  the  stronger  effect 
on  both  the  flow  and  the  heat  transfer  of  the  cavity.  Their  results  at  finite  Marangoni 
number  in  a  finite  pool  suggest  that  experimentally  observed  instabilities  may  be 
associated  with  the  rapid  turning  flows  and  high  vorticity  in  the  cold  wall  region. 
Another  important  finding  in  this  paper  was  the  determination  of  the  scalings  for  the 
thermocapillary  boundary-layer  thickness  on  the  free  surface  and  on  the  rigid  wall. 
Their  results  confirmed  those  of  Ostrach  [Ref.  5,  24]. 

Ohring  and  Lugt  [Ref.  10]  studied  the  onset  of  Marangoni  convection  in  a  float 
zone  of  liquid  silicon  from  a  state  at  rest  in  the  absence  of  gravity,  with  high  Marangoni 
numbers.  They  found  the  existence  of  a  critical  Marangoni  number  at  which  the 
steady  flow  field  becomes  unstable  and  changes  to  an  axisymmetric  oscillatory  field. 
Also,  they  concluded  that  a  deformable  free  surface,  even  a  very  small  one  on  the 
order  of  10~4  cm,  is  necessary  for  the  onset  of  instability.  The  velocity  at  the  free 
surface  is  extremely  high.  The  temperature  field  oscillates  and  causes  uneven  heat 
transfer  at  the  walls,  which  is  not  desirable  for  silicon  crystal  growth. 

Chan,  Mazumder,  and  Chen  [Ref.  18]  developed  a  three-dimensional  model  of 
the  fluid  flow  and  heat  transfer  of  a  laser  melted  pool.  Using  a  perturbation  solution, 
they  found  that  the  presence  of  the  thermocapillary  convection  causes  the  physics  of 
the  problem  to  change  from  conduction  to  convection  dominated.  The  pool  geometry 
then  changes  dramatically,  resulting  in  up  to  a  150  percent  increase  in  the  aspect  ratio 
as  compared  to  the  pure  conduction  case.  In  another  paper  [Ref.  12],  they  presented 
solutions  of  thermocapillary  convection  in  the  central  region  of  a  nonuniformly  heated 
surface  for  the  asymptotic  limits  of  very  high  and  very  small  Prandtl  numbers.  The 
model  is  valid  when  the  viscous  and  thermal  boundary  layers  are  small  compared  to 
the  depth  and  width  of  the  melt  pool.  In  addition,  the  intensity  of  the  thermocapil¬ 
lary  convection  and  the  boundary-layer  thicknesses  in  the  stagnation  region  depend 
primarily  on  the  curvature  of  the  heat  flux  distribution.  In  a  comprehensive  report, 
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Mazumder,  Chen,  Chan,  and  Zehr  [Ref.  7]  developed  a  three-dimensional  perturba¬ 
tion  model  for  surface  tension  driven  flow  with  a  flat  surface  which  predicted  velocity 
field,  temperature  field,  effect  of  trace  elements  and  a  self-consistent  prediction  of  the 
pool  shape  and  cooling  rate.  The  three-dimensional  problem  is  represented  by  two 
sets  of  two-dimensional  governing  equations.  It  is  found  that  a  particle  recirculates 
many,  many  times  in  the  molten  pool  before  it  resolidifies,  showing  that  the  solute 
can  be  well  mixed.  Regarding  free  surface  deformation,  thermocapillarity  drives  the 
surface  fluid  radially  outward  at  extremely  high  velocities.  These  high  velocities  dis¬ 
place  more  mass  from  the  surface  region  than  can  be  replaced  by  the  recirculating 
flow,  which  thus  causes  a  depression.  The  displaced  mass  builds  up  at  a  solid/liquid 
interface,  which  causes  the  surface  to  bulge  upward  where  the  liquid  turns  downward 
into  the  molten  pool.  Experiments  were  performed  to  examine  the  validity  of  the 
theoretical  models,  and  both  theoretical  models  and  experiments  predict  the  same 
trend  in  aspect  ratio:  the  aspect  ratio  increases  with  laser  power. 

Masud,  Kamotani,  and  Ostrach  [Ref.  19]  experimentally  investigated  high 
Prandtl  number  flow  in  a  half-zone  configuration.  Specifically,  they  studied  the  ef¬ 
fects  of  buoyancy  and  column  shape  on  the  onset  conditions  of  oscillations.  There 
is  a  minimum  critical  Marangoni  number,  below  which  no  oscillations  appear.  They 
speculated  that  both  convection  and  a  dimensionless  surface  deformation  parameter 
must  be  taken  into  account  together  to  explain  the  oscillation  mechanism.  Kamotani 
and  Ostrach  [Ref.  20]  then  conducted  a  theoretical  analysis  of  the  half-zone  con¬ 
figuration.  They  concluded  that  free  surface  deformation  plays  an  important  role 
in  oscillatory  thermocapillary  flow  in  high  Prandtl  number  fluids.  The  deformation 
induces  a  change  in  the  surface  flow  and  alters  the  driving  force  in  the  hot  corner. 
This  then  triggers  oscillation  cycles  in  which  the  flow  at  the  surface  alternates  be¬ 
tween  fast  and  slow.  Finally,  they  derive  a  surface  deformation  parameter  by  scaling 
analysis,  alluded  to  in  their  earlier  paper  [Ref.  19].  The  scaling  for  the  core  velocity 
found  by  Cowley  and  Davis  was  confirmed  for  several  high  Prandtl  number  flows  in 
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the  Kamotani  and  Ostrach  paper.  Kuhlmann  [Ref.  3]  points  out  that  this  behavior 
suggests  that  for  high  Prandtl  numbers  and  Reynolds  numbers  that  axe  not  too  large, 
the  hot  corner  can  determine  the  scalings  of  the  surface  temperature  and  the  core 
velocity. 

As  stated  earlier,  the  Marangoni  effect  is  caused  by  a  jump  in  shear  stress 
across  the  interface  which  balances  the  surface-tension  gradient  along  the  interface. 
Scriven  and  Sternling  [Ref.  13]  give  a  short  but  interesting  historical  review  of 
this  phenomenon,  including  causes  and  resistance  to  surface  movements.  Batischev, 
Kuznetsov,  and  Pukhnachov  [Ref.  21]  investigated  the  Marangoni  effect,  dividing 
the  cases  involving  the  thermocapillary  effect  into  classes  for  which  the  motion  shows 
similarity,  but  they  mainly  considered  the  case  where  the  Reynolds  number  is  high 
and  the  Prandtl  number  is  of  order  unity.  They  define  the  Marangoni  boundary  layer 
as  a  surface  layer  where  intensive  convection  quickly  decays  in  depth.  The  fluid  in 
the  Marangoni  boundary  layer  differs  in  a  qualitative  sense  from  that  in  the  Prandtl 
classic  boundary  layer,  which  is  initiated  by  an  external  fluid  flow  due  to  the  no-slip 
condition  at  a  rigid  wall.  The  Marangoni  boundary  layer  is  initiated  by  forces  acting 
on  a  free  surface  and  subsequently  initiates  the  fluid  flow  inside  the  entire  volume  [Ref. 
21].  Pukhnachov  [Ref.  22]  systemized  the  similarity  solutions  of  the  Navier-Stokes 
equations  in  another  paper,  distinguishing  and  justifying  boundary-layer  asymptotics 
when  the  viscosity  goes  to  zero,  describing  the  construction  of  nonstationary  bound¬ 
ary  layers  on  a  liquid  free  surface  and  providing  descriptions  of  asymptotic  solutions 
in  several  cases. 

Napolitano  [Ref.  23]  also  studied  Marangoni  boundary  layers  between  two 
immiscible  fluids  in  the  plane,  restricting  analysis  to  the  case  with  a  Prandtl  number 
of  order  unity.  Three  necessary  conditions  for  the  existence  of  a  Marangoni  boundary 
layer  were  given.  First,  the  region  must  be  sufficiently  far  away  from  solid  boundaries; 
second,  surface  driving  forces  must  be  of  the  same  order  of  magnitude  as  viscous  forces; 
and  third,  the  thickness  of  the  region  must  be  much  smaller  than  the  thickness  of  the 
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interface  between  the  two  immiscible  fluids.  The  influence  of  the  flow  field  on  the 
dynamic  interface’s  shape  was  also  calculated  in  terms  of  characteristic  length  ratios. 
Napolitano  [Ref.  23]  and  Ostrach  [Ref.  24]  separately  derived  the  scaling  law  of  low 
Prandtl  number  inertial  flow  in  a  free  surface  layer. 

Strani,  Piva,  and  Graziani  [Ref.  25]  also  studied  thermocapillary  convection 
in  a  rectangular  cavity.  They  found  an  asymptotic  solution  in  the  limiting  case  where 
the  aspect  ratio  goes  to  zero.  For  increasing  values  of  the  aspect  ratio,  the  motion 
is  confined  in  a  region  near  the  free  surface.  As  with  other  investigators,  Strani, 
Piva,  and  Graziani  concluded  that  any  surface  deformation  seems  to  have  a  negligible 
influence  on  the  qualitative  aspects  of  the  flow  field.  The  solution  method  employed, 
although  independently  derived,  turns  out  to  be  a  special  case  of  the  method  proposed 
by  Sen  and  Davis  [Ref.  11].  Regarding  the  cold  comer,  they  concluded  that  the  larger 
surface  temperature  gradient  near  the  cold  wall  develops  a  local  increase  in  the  driving 
force.  This  leads  to  a  complex  stagnation  flow  field,  with  significant  acceleration  near 
the  wall,  positive  until  the  maximum  velocity  is  reached,  then  negative  to  meet  the 
boundary  condition  of  zero  velocity.  A  boundary-layer-type  flow  grows  downward 
along  the  cold  wall,  starting  from  the  stagnation  point  at  the  surface. 

Canright  [Ref.  1]  introduced  study  of  the  dynamics  restricted  solely  to  the  cold 
corner  region,  where  the  flow  along  the  free  surface  toward  the  cold  wall  compresses 
the  thermal  gradient,  thereby  enhancing  the  flow  in  a  sort  of  positive  feedback.  This 
is  known  as  the  cold  corner  problem.  This  feedback  results  in  small  local  length 
scales  and  high  velocities  near  the  corner.  Numerical  results  are  consistent  with 
scaling  analysis.  In  particular,  for  frilly  convective  flow,  the  comer  behavior  is  lo¬ 
cally  determined.  Scaling  limits  for  the  Reynolds  and  Marangoni  numbers  are  given 
for  four  different  regimes:  conductive- viscous,  convective-viscous,  conductive-inertial, 
and  convective-inertial.  Canright  concluded  that  increasing  the  global  Marangoni 
number  decreases  the  local  length  scale  to  give  an  effective  local  Marangoni  number 
of  unity  [Ref.  1].  In  contrast  to  the  hot  corner  problem  posed  by  Cowley  and  Davis, 
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Canright  found  that  no  thermal  boundary  layer  forms  when  the  Prandtl  number  is 
large.  This  is  due  to  the  surface  forcing  being  limited  to  a  relatively  concentrated 
region  in  the  cold  corner,  while  for  the  hot  corner,  the  forcing  is  distributed  over  a 
broad  region,  since  the  thermal  variations  in  the  horizontal  direction  are  extended  by 
convection.  Canright  also  validated  his  scales  by  comparing  his  results  with  those  of 
Zebib,  Homsy,  and  Meiburg  [Ref.  17],  who  showed  that  the  temperature  gradient  on 
the  free  surface  can  vary  considerably  within  a  small  distance  from  either  of  the  hot 
and  cold  walls.  Canright’s  scaling  of  the  boundary  layers  were  also  published  in  a 
treatment  of  the  cold  corner  by  Kuhlmann  [Ref.  3]. 

Huber  [Ref.  26]  extended  the  study  of  the  cold  corner  problem,  examining 
the  problem  numerically  for  different  low  Marangoni  numbers  using  a  Green’s  func¬ 
tion  flow  representation  for  the  viscous  case,  in  the  limit,  as  the  Reynolds  number 
approaches  zero.  One  advantage  of  the  Green’s  function  is  that  the  flow  can  still  be 
represented  over  the  entire  quarter  plane,  so  that  there  is  no  artificial  recirculation 
due  to  imposed  artificial  boundaries.  The  Green’s  function  approach  does  not  require 
flow  boundary  conditions  at  the  computational  domain’s  boundaries.  The  flow  is 
assumed  isothermal  across  the  boundary  and  is  recirculating,  decaying  with  distance. 

Canright’s  work  in  1994  concentrated  on  the  Marangoni  number  dependence 
of  the  cold  corner  problem.  A  representative  velocity  vector  field  (using  Canright’s 
numerical  data)  is  shown  in  Figure  1.  This  work  reconsiders  the  problem,  refining 
the  Prandtl  number  dependence  of  the  variables.  All  the  graphs  shown  in  [Ref.  1] 
are  correct,  but  we  will  now  show  that  the  theoretical  scales  in  the  cold  corner  for 
the  convective  inertial  case  are  different  than  those  previously  published,  in  terms 
of  the  Prandtl  number  dependence.  Canright’s  numerical  results  reflect  a  boundary- 
layer  structure  that  we  assume  throughout  this  work  and  use  as  a  guiding  principle. 
The  task  is  to  solve  for  the  flow  and/or  heat  in  each  region  and  match  it  to  that  in 
neighboring  regions  in  order  to  predict  future  flow  conditions.  The  scalings  derived 
herein  are  different  than  those  published  by  Canright;  nonetheless,  they  agree  better 
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Distance  From  Free  Surfa< 


Figure  1.  Velocity  Vector  Field 


with  his  numerical  data. 
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III.  PROBLEM  STATEMENT 


Most  of  the  fundamental  ideas  of  science  are  essentially  simple,  and  may,  as 

A  RULE,  BE  EXPRESSED  IN  LANGUAGE  COMPREHENSIBLE  TO  EVERYONE.  — Albert  Einstein 


A.  MATHEMATICAL  FORMULATION  OF  THE  PROB¬ 
LEM 

A  pool  of  incompressible  Newtonian  fluid  is  bounded  on  the  left  by  a  vertical 
solid  wall,  maintained  at  a  cold  temperature  Tc,  while  the  undisturbed  fluid  far  from 
the  cold  corner  is  at  the  hot  ambient  temperature  Th  of  the  core  flow  (see  Figure  2 
below). 

AIR 


Figure  2.  Problem  Formulation 


Above  the  horizontal  free  surface  of  the  liquid  is  an  inviscid,  nonconducting 
gas.  Surface  tension  is  assumed  strong  enough  to  keep  the  free  surface  flat  (small 
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Capillary  number),  but  with  surface  tension  variations  due  to  a  linear  dependence  on 
temperature.  The  resulting  flow  is  assumed  to  be  two-dimensional  and  steady. 

The  equations  governing  the  thermocapillary  convection  in  the  cold  corner  are 
the  conservation  of  mass,  the  conservation  of  momentum,  and  the  conservation  of 
energy  equations: 

V  •  u  =  0  (III.l) 

pu-Vu  =  -Vp  +  pV2u  (III. 2) 

pcpU-  VT  =  kV2T  (III.3) 

with  the  boundary  conditions  that  at  y  =  0,  Ty  =  0,  v  —  0,  and  puy  =  ^Tx.  At 
x  =  0,  T  —  Tc,  and  u  =  v  =  0.  As  (x,  y)  — ►  oo,  T  — >  Th,  and  u,  v  — >  0. 

In  these  equations,  u  is  the  velocity  vector  with  components  u  and  v  in  the 
x  (horizontally  rightward)  and  y  (vertically  downward)  directions,  p  is  the  pressure, 
T  is  the  temperature,  p  is  the  density,  p,  is  the  viscosity,  Cp  is  the  specific  heat,  k 
is  the  thermal  conductivity,  and  the  surface  tension  is  assumed  to  be  of  the  form 
<7  =  <rc  —  7  (T  —  Tc),  with  7  a  positive  constant.  The  boundary  conditions  specify 
that  the  cold  wall  is  isothermal  with  fluid  obeying  the  no-slip  condition,  and  the  flat 
free  surface  is  thermally  insulated,  with  thermocapillary  forcing. 

To  nondimensionalize  the  equations,  a  heat  flux  scale  of  Q  and  a  temperature 
scale  (relative  to  the  cold  temperature)  of  AT  =  Th~Tc  axe  introduced.  Thermal  con¬ 
duction  gives  the  length  scale  of  d  =  Q/k  AT,  the  thermocapillary  coupling  gives  the 
velocity  scale  u  =  7 A T / p,  and  the  viscous  pressure  scale  is  p  =  pu/d  =  k-y  (AT)2/ Q. 
In  summary,  the  nondimensional  scales  are: 


X  = 

(Q/k  AT)  x' 

(III.4) 

u  = 

(7  AT /p)  u' 

(III.5) 

p  = 

(k1(ATf/Q)p> 

(HI-6) 

T  = 

(T-Th)/(TC-Th) 

(III.7) 
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which  implies  T  —  T^  —  (AT)  T1,  where  the  variables  with  primes  axe  dimensionless. 
Begin  with  the  conservation  of  energy  equation,  which  can  be  rewritten  as 

u-VT  =  kV2T  (III.8) 

where  k  is  the  thermal  diffusivity.  Substituting  the  dimensionless  scales  into  this  form 
of  the  energy  equation  yields 

^[u'.VT']  =  V2r  (III.9) 

Dropping  the  primes  and  defining  a  dimensionless  parameter  M,  the  Marangoni  num¬ 
ber,  asM  =  ^|,  the  nondimensional  conservation  of  energy  equation  becomes 

M(u-VT)  =  V2T.  (III.  10) 

Next,  the  conservation  of  momentum  equation  may  be  rewritten  as: 

u  ■  Vu  =  — ^  Vp  +  v  V2u  (III.  11) 

where  v  is  the  kinematic  viscosity.  Substituting  the  dimensionless  scales  into  the 
momentum  equation  yields 

p  •  w]  =  -Vp'  +  VV.  (III.  12) 

Dropping  the  primes  and  defining  a  dimensionless  parameter  R,  the  Reynolds  number, 
as  R=  the  nondimensional  conservation  of  momentum  equation  becomes 

R  (u  •  Vu)  =  -Vp  +  V2u.  (III.  13) 

Situations  in  which  the  Reynolds  number  is  small  are  called  slow  viscous  flows,  because 
the  viscous  forces  arising  from  shearing  motions  of  the  fluid  predominate  over  inertial 
forces  associated  with  the  acceleration  or  deceleration  of  the  fluid  particles.  The 
thermocapillary  boundary  condition  is  also  nondimensionalized  by  substituting  (III.4) 
and  (III.5)  into  puy  =  7 Tx ,  obtaining  (dimensionless)  uy  =  Tx.  In  addition,  the  mass 
conservation  equation,  V  •  u  =  0,  remains  unchanged  after  nondimensionalization. 
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The  ratio  of  the  Marangoni  number  and  the  Reynolds  number  gives  the  nondi- 
mensional  Prandtl  number 

P  =  v/k  =  M/R.  (III.  14) 

We  want  to  analyze  the  system  for  large  values  of  the  Marangoni  number,  M,  which 
measures  the  importance  of  thermal  convection  relative  to  thermal  diffusion.  For  large 
enough  M,  convection  becomes  important  and  the  strong  surface  flow  towards  the 
cold  wall  compresses  the  thermal  gradient  along  the  surface,  which  in  turn  strengthens 
the  local  driving  force  for  the  flow.  Also,  since  R  =  M/P ,  the  effect  of  making  R  1 
is  to  make  the  inertia  force  much  greater  than  the  viscous  force,  so  that  inertia  and 
pressure  forces  are  dominant,  except  along  the  rigid  wall  and  free  surface,  where 
viscous  boundary  layers  are  formed. 

Near  the  surface,  vorticity  must  be  taken  into  account.  Vorticity  is  the  curl 
of  the  velocity,  Q  =  V  x  u,  which  (in  two-dimensional  flow)  reduces  to  the  one 
component  along  the  axis  normal  to  the  x-y  plane,  yielding  u  =  vx  —  uy.  Thus, 
eliminating  pressure  by  taking  the  curl  of  the  momentum  equation,  the  following 
equation  is  obtained: 

u  •  Vu)  =  is  V2u>  (III.  15) 

This  equation  is  referred  to  as  the  vorticity  transport,  or  vorticity  transfer, 
equation.  It  states  that  the  substantive  variation  of  vorticity,  which  consists  of  the 
local  and  convective  terms,  is  equal  to  the  rate  of  diffusion  of  vorticity  through  friction. 
In  our  non-dimensionalized  form,  the  vorticity  equation  becomes  R  (u  ■  Vu;)  =  V2 u. 
At  the  surface,  the  thermocapillary  condition  gives  uj  =  Tx,  which  comes  from  u  = 
vx  —  uy,  and  the  fact  that  vx  =  0  at  the  surface. 

B.  ASSUMPTIONS  AND  SCALING  ANALYSIS 

A  few  assumptions  concerning  the  scaling  analysis  are  made.  First,  a  high 
Reynolds  number,  R,  implies  that  inertial  forces  are  dominant,  outside  the  viscous 
boundary  layers.  Second,  a  high  Marangoni  Number,  M,  implies  that  thermal  con- 
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vection  is  important.  Third,  a  low  Prandtl  Number,  P,  implies  that  there  exists  a 
low  ratio  of  viscous  to  thermal  diffusion  and  the  region  is  inertial. 


AIR  Ty  =  o 


Figure  3.  Scaling  of  the  Domain 


Recall  that  we  are  assuming  the  structure  based  on  Canright’s  numerical  data. 
Figure  3  shows  a  schematic  of  the  corner  with  viscous  boundary  layers  along  the  free 
surface  and  rigid  wall.  The  depth  and  width  of  the  pool  are  large  compared  to  all 
local  length  scales,  so  the  pool  appears  semi-infinite  both  horizontally  and  vertically. 

To  distinguish  the  various  notations,  a  hat  Q  will  denote  the  original  nondi- 
mensionalized  variables.  Capital  letters  will  denote  the  core  scaling,  and  lowercase 
letters  will  denote  boundary-layer  variables  that  do  not  have  the  same  scaling  as  the 
core  variables. 

Let  the  core  velocity  scale  be  denoted  by  U.  Let  the  characteristic  velocity 
scale  along  the  free  surface  be  u.  We  assume  that  surface  thermal  variations  will 
be  confined  to  a  region  with  a  horizontal  characteristic  length  scale  of  l.  We  also 
assume  that  vorticity  will  be  confined  to  thin  viscous  boundary  layers.  The  flow  in 
the  core  region  is  dominated  by  inertia.  (This  “defines”  the  core  region  as  the  area 
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with  negligible  vorticity).  The  solid  boundary  acts  as  a  source  of  vorticity,  which  is 
then  diffused  away  by  viscosity  and  converted  downstream  with  the  fluid. 

Along  the  surface,  the  boundary  layer  results  from  thermocapillary  stresses 
and  will  have  thickness  8.  Along  the  wall,  the  boundary  layer  results  from  the  no-slip 
condition  and  will  have  thickness  A.  Thus,  8  is  the  vertical  length  scale  of  velocity 
shear  along  the  free  surface,  and  A  is  the  horizontal  length  scale  of  velocity  shear 
along  the  rigid  wall. 

In  case  of  large  M,  convection  is  important  and  the  rapid  surface  flow  into 
the  cold  corner  compresses  the  thermal  gradient  along  the  surface,  reducing  l.  The 
length  l  will  be  reduced  to  the  point  that  thermal  diffusion  away  from  the  rigid  wall 
will  balance  the  strong  convection  toward  the  wall. 

In  case  of  large  R,  inertia  is  dominant  and  there  will  be  two  viscous  boundary 
layers  of  thickness  8  and  A  (~  8)  along  the  surface  and  rigid  wall,  respectively.  The 
effect  of  increasing  the  Reynolds  number  to  values  large  compared  with  unity  is  to 
confine  the  vorticity  diffused  from  the  rigid  wall  to  a  layer  of  relatively  small  thickness. 
As  R  increases,  8,  the  thickness  of  the  boundary  layer,  decreases.  Convection  contains 
the  vorticity  inside  the  boundary  layer. 

Surface  tension  variations  are  due  to  an  assumed  linear  dependence  on  the 
temperature.  For  the  boundary  layer  along  the  free  surface,  the  thermocapillary 
stress  condition  at  the  surface  gives  uy  ~  Tx.  This  provides  the  relation  |  ~  or, 
solving  for  the  surface  velocity  in  the  horizontal  direction,  u  ~  f.  The  continuity 
condition  V  •  u  =  0  gives  the  relation  ux  ~  vy.  Therefore,  j  ~  or  v  ~  j,  which 
gives  v  ~  (j)2  as  the  vertical  component  of  the  velocity  near  the  surface. 

At  the  boundary-layer  edge/core  region  interface,  normal  velocities  must  match, 
so  v  ~  V,  and  continuity  requires  that  Vy  ~  Ux.  Thus  Vy  ~  ^  ~  Ux  ~  j,  which 
gives  the  relation  ^  ~  f  • 

Next,  a  check  of  the  dominant  balance  of  terms  is  required.  In  the  vorticity 
equation,  R  (jp  + Ignore  ^  as  too  small,  since  l  3>  8.  This  implies 


1  f  o 

that  R  J|  ~  jfi,  or  that  R~-  ~  1.  Therefore,  this  gives  an  expression  for  u  in  terms 
of  the  Reynolds  number,  u  ~  or  w  ~  jiM- 

Similarly,  in  the  energy  equation  applied  to  the  core  region,  MUTX  ~  Txx  (this 
dominant  balance  will  be  described  in  detail  later;  for  now,  assume  Ty  and  Tyy  scale 
comparably  to  Tx  and  Txx  in  the  core  flow  region).  Thus,  MUTX  ~  Mj  ~  Txx  ~ 
so  that  ~  -fa,  which  yields  another  relationship  for  the  surface  velocity,  u  ~ 

We  assume  that  the  wall  boundary  layer  is  passive,  fed  by  the  mass  flux  from 
the  surface  layer  (due  to  conservation  of  mass).  We  can  derive  its  scales  by  using  the 
mass  flux,  continuity  and  vorticity  equations.  Conservation  of  mass  implies  that  the 
velocity  out  of  the  corner  is  the  same  as  that  into  the  corner  (v  ~  u).  The  momentum 


equation  gives  R  (Uvx  +  Wy)  —  vxx  +  vyy,  and  continuity  gives  Ux  ~  vy.  This  means 
K  ~  l ,  or  v  ~  Substituting,  f  (U%  +  leaving 

A  ~  M~lPU~l  and  L  ~ 


To  obtain  scales  for  the  lengths  and  velocity,  use  the  three  relationships  for 
the  free  surface  horizontal  velocity  scaling:  Solving  for  the 

horizontal  boundary  layer  length  /,  l  ~  62M  ~  J;,  which  implies  that  <5  ~  Thus, 
l  ~  7^2,  «  ~  P,  and  U  ~  P2.  Further,  A  ~  ^  ~  6  and  L  ~  ^ 

To  summarize,  using  the  boundary  condition  at  the  free  surface  gives  a  starting 
point  for  the  scalings  and  dominant  balance.  It  is  important  to  obtain  relations  for 
the  thermal  length  scale  l,  the  viscous  thickness  6 ,  the  surface  velocity  u ,  and  the 
core  velocity  U  scales  in  the  cold  corner,  in  terms  of  the  dimensionless  parameters  M 
and  P.  These  scalings  are 


6  ~  M-1P-1, 
l  ~  M~lP~ 2, 
u  ~  P, 

U  ~  P2. 
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These  scalings  are  consistent  with  the  numerical  results  obtained  by  Canright  in  his 
1994  work  [Ref.  1]. 
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IV.  THE  COLD  CORNER  REGIONS 


It  is  not  knowledge,  but  the  act  of  learning,  not  possession,  but  the  act  of 

GETTING  THERE,  WHICH  GRANTS  THE  GREATEST  ENJOYMENT.  — Carl  Friedrich  Gauss 


A.  OVERVIEW 

The  flow  into  the  corner  is  assumed  to  be  convective  and  inertial,  due  to  the 
large  Marangoni  number  and  small  Prandtl  number  (M  »  1  ,P  <C  1),  governed  by 
the  system  of  equations 

R  (a- va)  =  V2£  (IV.  1) 

V-fi  =  0  (IV.2) 

M  (Z-VT)  =  V2f  (IV.3) 

The  cold  corner  is  divided  into  four  regions:  an  outer  core  region  away  from 

the  surface  and  cold  wall  where  the  flow  is  inviscid  and  relatively  simple;  an  inner 
viscous  boundary  layer  along  the  liquid  free  surface;  an  inner  viscous  boundary  layer 
along  the  wall/liquid  interface;  and  a  thermal  layer  overlapping  the  viscous  layers  and 
the  core  flow.  A  representative  velocity  vector  field  (using  Canright’s  numerical  data) 
was  shown  in  Figure  1  in  Chapter  II.  The  task  is  to  solve  for  the  flow  and/or  heat 
in  each  region  and  match  it  to  that  in  neighboring  regions  in  order  to  predict  future 
flow  conditions.  The  scaling  factors  are  6  ~  A  ~  M~lP~l,  l  Jj  M~lP~ 2,  u  ~  P, 
and  U  ~  P2.  In  the  core  flow  region,  let  the  horizontal  and  vertical  components  of 
the  velocity  be  designated  by  capital  letters.  The  governing  equations  in  each  region 
will  be  examined,  using  the  appropriate  nondimensionalized  scalings. 
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B.  VISCOUS  SURFACE  BOUNDARY  LAYER 

A  viscous  boundary  layer  forms  along  the  surface,  with  thickness  6  and  length 
l.  The  scalings  used  in  this  region  are  u  =  Pu,  v  =  P2V,  x  =  M~lP~2X,  y  = 
y,  T  =  T  +  •  •  •,  and  p  =  M  P3p.  The  pressure  scales  with  the  dominant 
term  in  the  core  flow.  The  momentum  equation  in  the  x-direction  can  be  written  as 

uux  +  Vuy  =  —P2px  +  P 2  uxx  +  Uyy,  (IV.4) 

where  to  leading  order  we  can  neglect  the  terms  of  order  P 2 .  The  boundary  conditions 
at  y  =  0  give  uy  =  Tx  and  V  =  0.  The  continuity  equation  becomes 

ux  +  Vy  =  0.  (IV.5) 

Prom  the  momentum  equation  in  the  y-direction,  the  pressure  change  in  the  y- 
direction  is  negligible  through  the  boundary  layer,  indicating  that  the  pressure  is 
effectively  constant  throughout  the  layer  for  each  X  value.  Mathematically,  px  py- 
It  is  assumed  then  that  the  pressure  distribution  in  the  boundary  layer  is  imposed  by 
the  core  flow  pressure  gradient  just  outside  the  boundary  layer.  In  Ludwig  Prandtl’s 
words,  “the  pressure  distribution  will  be  impressed  on  the  transition  layer  by  the  free 
fluid”  [Ref.  27].  Also,  the  velocity  in  the  free  surface  boundary  layer  is  continuous 
with  the  core  velocity  U. 

The  energy  equation  in  the  surface  viscous  boundary  layer  is  derived  using  a 
boundary-layer  expansion  for  T  and  will  be  given  later. 

C.  VISCOUS  WALL  BOUNDARY  LAYER 

A  viscous  boundary  layer  forms  along  the  wall,  with  thickness  A  and  length 
L.  This  viscous  boundary  layer  is  different  from  the  surface  viscous  boundary  layer 
in  that  there  is  a  no-slip  condition  at  the  wall.  The  scalings  used  in  this  region  are 
u  =  P2U,v  =  Pv,x  —  M~1P~1  x,  y  =  M~lP~ 2  Y,  T  =  T  +  •  •  •,  and  p  =  M P3 p. 
The  momentum  equation  in  the  y-direction  can  be  written  as 

U  Vx  +  V  Vy  —  — P2  py  +  Vxx  +  P2  Vyy.  (IV.6) 
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where  again  terms  of  order  P2  can  be  neglected  to  leading  order.  The  no-slip  condition 
at  the  rigid  wall  gives  u  =  v  =  0atx  =  0.  Similarly,  the  continuity  equation  becomes 

Ux  +  vY  =  0.  (IV.7) 

The  energy  equation  for  the  viscous  boundary  layer  on  the  wall  is 

UTx  +  vTy  =  p-1  Txx  +  P  Tyy  (IV.8) 

with  the  boundary  condition  at  x  =  0  of 

T  =  -1  (IV.9) 

The  P  Tyy  term  is  neglected  as  small.  The  Txx  term  must  therefore  be  small  enough 
to  balance  the  equation.  Thus,  there  is  no  need  to  solve  the  energy  equation  in  this 
region,  as  the  temperature  is  basically  constant  across  the  wall’s  viscous  layer.  Also, 
due  to  the  conditions  of  high  Marangoni  number  coupled  with  low  Prandtl  number, 
the  larger  thermal  layer  encompasses  the  viscous  boundary  layer  at  the  wall. 

D.  CORE  FLOW  REGION 
1.  Flow  Equations 

The  core  flow  is  assumed  to  be  irrotational,  which  means  that  there  is  no 
vorticity,  and  also  is  incompressible,  which  means  that  density  is  not  affected  by 
changes  in  the  pressure.  When  the  velocity  U  is  solenoidal,  we  may  introduce  a 
stream  function  ip  which  allows  us  to  replace  the  horizontal  and  vertical  components 
of  velocity  by  a  single  function.  Define  ip  by  u  =  and  v  =  — 1|.  These  defin¬ 
itions  satisfy  the  continuity  equation.  When  the  flow  is  also  irrotational,  then  the 
expression  of  zero  vorticity  gives  Laplace’s  equation,  V2ip  =  0.  Although  the  equa¬ 
tion  of  motion  is  nonlinear  in  U,  the  velocity  distribution  is  completely  determined 
by  a  linear  equation  derived  from  the  restrictive  condition  of  irrotationality  and  the 
mass-conservation  equation.  The  nonlinear  equation  of  motion  is  needed  only  for  the 
calculation  of  the  pressure  after  the  velocity  distribution  has  been  determined.  The 
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boundary  conditions  for  the  core  flow  must  match  the  normal  velocities  in  each  of 
the  boundary  layers.  Scale  the  velocity  in  this  region  by  P2.  Scale  the  horizontal  and 
vertical  lengths  by  M~lP~2.  Scale  the  pressure  by  MP3.  The  momentum  equation 
in  the  ^-direction  is 

U  Ux  +  VUy  =  —px  +  P  ( Uxx  +  Uyy)-  (IV.  10) 

The  momentum  equation  in  the  y-direction  is 

UVx  +  VVy  =  -pY  +  P(Vxx  +  Vyy).  (IV.ll) 

With  M  1,  and  P  «  1,  these  momentum  equations  reduce  to  U  Ux  +  V  Uy  =  ~Px 
in  the  ^-direction  and  UVx  +  VVy  =  —py  in  the  y-direction.  The  momentum 
equations  could  also  be  written  in  terms  of  the  stream  function,  ip.  For  example,  the 
momentum  equation  in  the  x-direction  becomes 

1pY  IpXY  ~  4>X  'pYY  =  —  PX  +  P('lPxXY  +  'tpYYY )  (IV.  12) 


Laplace’s  equation  can  be  written  as 


2 d2ip  ,  d2ip 


Vzip  = 


+ 


=  0 


(IV.  13) 


dx 2  dY2 

In  a  two-dimensional  flow,  the  lines  of  constant  ip  are  known  as  the  streamlines,  and 
the  difference  between  the  numerical  values  of  two  streamlines  is  equal  to  the  flow  rate 
between  the  streamlines.  An  equivalent  description  of  two-dimensional  incompressible 
irrotational  flow  is  in  terms  of  a  velocity  potential  (p ,  such  that  u  =  X(p.  The  stream 
function  ip  is  harmonic,  as  is  0,  and  together  they  satisfy  the  Cauchy-Riemann  equa¬ 
tions  as  conjugate  harmonic  functions,  which  means  that  lines  of  constant  ip  and  (p 
are  orthogonal.  The  fact  that  ip  and  (p  are  harmonic  and  satisfy  the  Cauchy-Riemann 
equations  is  necessary  and  sufficient  for  the  definition  of  a  complex  potential  F,  de¬ 
fined  by  F  =  <p  +  iip,  where  z  =  x  +  iy  and  the  stream  function  and  velocity  potential 
are  both  functions  of  x  and  y.  F  is  analytic,  and  note  that  F'(z )  =  (px(x,  y)  +iipx(x,  y), 


or  F'(z)  =  <px(x,  y)  -  i<py(c r,  y).  The  velocity  thus  becomes  u  =  F'(z),  and  the  speed, 
or  magnitude  of  the  velocity,  is  found  by  |u|  =  \F'{z)\. 

Since  the  equation  V2ip  =  0  is  linear,  we  may  superpose  solutions  for  different 
flows  and  add  directly  the  values  of  ip  at  every  point  in  the  plane  to  obtain  the  new 
values  of  ip,  which  represent  physically  the  direct  superposition  of  the  various  flows. 

2.  Dominant  Balance 

Naive  APPLICATION  OF  THE  METHOD  OF  DOMINANT  BALANCE  IS  BOOTLESS.  —  Carl  M. 

Bender  and  Steven  A.  Orszag 


How  does  the  temperature  T  change  across  the  boundary  layer?  Let  U  be  the 
core  velocity  and  T  be  the  core  temperature,  and  u  and  t  are  the  additional  values 
of  the  velocity  and  temperature  inside  the  surface  viscous  boundary  layer,  so  u  and 


t  are  significant  only  inside  the  layer.  The  horizontal  derivative  of  the  temperature, 
Tx,  in  the  core  flow  region  is  similar  to  Tx  at  the  surface. 


Let  T  =  T(X,Y )  +  rt(X,y),  where  r  is  an  unknown  scale.  Further,  let 
u  =  P2U(X,Y )  +  Pu(X,y),  using  the  known  scales:  x  =  ,  y  =  . 

To  determine  the  size  of  r,  substitute  the  expressions  for  T  and  u  into  the  energy 
equation.  This  yields 


M  [(. P2U  +  Pu){T  +  rt)£  +  P2V(T  +  rt) t]  =  (T  +  rt)±i  +  (T  +  rt)yy  (IV.14) 
or,  in  expanded  form, 


M2PaUTx  +  M2P\Tx  +  M2PiUrtx  +  M2Pzurtx  +  M2PAVTy  +  M2P3Vrty 

=  M2P4Txx  +  M2P4rtxx  +  M2PaTyy  +  M2P\tyy 

Divide  this  equation  by  M2PA  and  subtract  the  core  flow  region  equation,  UTX  + 
VTy  =  Txx  +  Tyy,  which  leaves 

-puTx  +  rUtx  +  —rutx  +  ^prVty  =  rtxx  +  rtyy.  (IV.15) 
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What  is  the  dominant  balance?  That  is,  how  must  the  scale  r  be  chosen  to  balance 
the  additional  heat  convection  of  the  first  term  in  (IV.  15)? 

Suppose  p uTx  ~  rUtX-  Then  r  ~  V,  which  implies  the  next  term  -p rutx  is 
0(p 2)-  Since  that  is  much  larger,  the  assumed  balance  cannot  be  dominant.  Thus, 
the  assumption  is  not  consistent. 

Suppose  puTx  ~  p rutx ■  Then  r  is  0(1),  and  the  term  p^rtyy  is  larger,  which 
means  the  assumed  balance  is  not  dominant  and  the  assumption  is  not  consistent. 

Suppose  puTx  ~  prVty.  This  is  equivalent  to  the  previous  case,  so  it  is  not 
consistent. 

Suppose  -puTx  ~  rtxx-  Again,  this  would  mean  that  r  ~  p,  as  in  the  first 
case,  which  is  not  consistent. 

Suppose  puTx  ~  j&Ttyy.  Then  r  ~  P  and  all  the  other  terms  are  small  in 
comparison  to  the  dominant  terms.  This  is  the  only  consistent  case.  Hence  r  —  P. 

Thus,  the  change  in  temperature  T  across  the  boundary  layer  is  Pt,  i.e.  O(P), 
so  it  is  negligible.  This  implies  that  the  temperature  gradient  Tx  at  the  top  of  the 
core  is  equivalent  to  that  right  on  the  surface,  which  drives  the  boundary-layer  flow. 
But  the  vertical  temperature  gradient  Tg  =  MP 2  [TV  +  ty],  so  the  surface  boundary 
condition  Ty  —  0  implies  TV  =  ~ty  at  y  =  0.  The  heat  equation  for  the  surface 
boundary  layer  is 

UTX  =  tyy  (IV.  16) 

In  order  to  determine  how  surface  layer  convection  modifies  heat  flux  into  the  core, 
use  the  boundary  condition  Ty  =  0  and  integrate  the  left-hand-side  of  the  surface 
heat  equation  across  the  surface  boundary  layer.  As  y  — >■  00,  ty  — >  0  and  the  heat 
flux  condition  in  the  core  becomes 

Ty(X,  0)=  f°°  uTx  dy.  (IV.17) 

Jo 
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Since  /0°°  u  dy  is  just  the  mass  flux  into  the  corner  (or  the  expression  for  the  stream 
function,  xp),  the  heat  flux  condition  becomes 


TY(X,0)=Txxfi(X,y->oo). 


(IV.  18) 
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V.  APPROXIMATE  METHODS 


An  idea  which  can  be  used  once  is  a  trick.  If  it  can  be  used  more  than  once 

IT  BECOMES  A  METHOD.  — George  Polya  and  Gabor  Szego 


The  goal  of  this  chapter  is  to  match  the  velocity  in  the  surface  viscous  boundary 
layer  to  that  in  the  wall  viscous  boundary  layer  and  to  that  in  the  core  flow  region. 
The  temperature  can  then  be  found  using  the  core  velocity.  For  the  viscous  boundary 
layers,  several  approximate  methods  exist  based  on  integrated  forms  of  the  boundary- 
layer  equations.  These  are  integral  methods  which  do  not  attempt  to  satisfy  the 
governing  equations  of  the  boundary  layer  for  each  streamline.  Instead,  the  equations 
are  satisfied  on  average  over  the  thickness  of  the  boundary  layer.  First,  start  with 
the  viscous  boundary  layer  along  the  free  surface,  since  the  surface  thermal  gradient 
drives  the  flow  into  the  comer. 

A.  VISCOUS  SURFACE  BOUNDARY  LAYER 
1.  Background 

The  first  approximate  method  considered  is  that  developed  by  Pohlhausen  in 
1921  [Ref.  28].  It  is  based  on  integrating  the  momentum  equation  across  a  boundary 
layer.  The  integrated  momentum  equation  gives  an  ordinary  differential  equation 
for  the  thickness  of  the  boundary  layer,  S(x),  provided  that  a  suitable  form  for  the 
velocity  profile  is  assumed.  Rosenhead  [Ref.  29]  and  Schlichting  [Ref.  30]  offer  two 
excellent  overviews  of  this  method,  as  well  as  descriptions  of  updated  methods  in  more 
modern  forms.  The  key  to  the  approximate  methods  is  that  they  allow  calculation 
of  the  momentum  thickness,  the  displacement  thickness,  and  the  shearing  stress  at 
a  wall.  These  will  all  be  defined  later.  Most  of  the  approximate  methods  for  two- 
dimensional  flows  in  the  literature  deal  with  flow  past  a  rigid  wall.  No  such  method 
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for  flow  along  a  free  surface  (where  the  no-slip  condition  of  the  rigid  wall  does  not 
apply)  has  yet  been  discovered  in  the  literature. 

What  is  the  critical  boundary-layer  thickness,  8(x)?  Pohlhausen’s  method  de¬ 
termines  the  boundary-layer  thickness  by  first  writing  the  momentum  equation  in 
terms  of  the  shear  stress,  the  displacement  thickness,  and  the  momentum  thickness 
of  the  boundary  layer.  The  displacement  thickness  8\  has  a  simple  physical  interpre¬ 
tation.  U6\  is  the  decrease  of  the  mass  flux  across  a  normal  to  the  surface,  due  to  the 
boundary  layer.  The  streamlines  of  the  outer  flow  are  thus  displaced  away  from  the 
surface  through  a  distance  6\  [Ref.  29].  It  is  the  distance  that  a  wall  would  have  to 
be  displaced  outward  into  the  free  stream  in  order  not  to  change  the  flow  field  if  the 
fluid  were  completely  inviscid  and  there  were  no  boundary  layer.  The  fluid  becomes 
slowed  down  by  viscosity  near  the  surface  and  is  thus  forced  outward,  so  that  the 
effective  surface  presented  to  the  oncoming  stream  is  increased  (or  thickened)  by  the 
boundary  layer’s  displacement  thickness  [Ref.  32].  The  momentum  thickness  has  the 
same  physical  significance  except  that  it  is  based  on  momentum  flux  instead  of  mass 
flux. 

Pohlhausen’s  approach  introduces  a  dimensionless  quantity  A,  a  shape  factor, 
which  can  be  interpreted  physically  as  the  ratio  of  pressure  forces  to  viscous  forces. 
According  to  Schlichting  [Ref.  30],  in  order  to  obtain  a  quantity  with  real  physical 
significance,  the  thickness  8  must  be  replaced  by  a  length  which  itself  has  physical 
significance,  such  as  the  momentum  thickness.  The  fundamental  idea  of  approximate 
methods  based  on  the  momentum  thickness  is  to  assume  that  the  dependence  of 
the  solution,  u(x,  y),  on  y  can  be  expressed  by  some  known  expression  of  y,  in  which 
coefficients  appear,  which  are  unknown  functions  of  x.  As  the  value  of  y  gets  large  (the 
distance  from  the  free  surface  increases),  the  boundary-layer  velocity  u  approaches 
the  core  velocity  U  outside  the  boundary  layer,  and  its  derivatives  tend  to  vanish. 
Generally,  only  as  many  boundary  conditions  for  y  =  0  are  taken  into  account  as  is 
necessary  to  determine  the  unknown  coefficients  in  the  velocity  profile.  This  will  be 
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done  later  in  this  section. 

In  view  of  this,  a  more  accurate  method  for  the  present  situation  is  that  of 
Timman  [Ref.  31],  whose  method  is  a  variant  of  Pohlhausen’s.  In  1949,  Tim  man 
developed  a  family  of  profiles  to  model  the  velocity  near  a  no-slip  surface  more  ad¬ 
equately  than  Pohlhausen.  Timman’s  method  seeks  to  find  an  exponential  form  for 
the  velocity  profile  u(x,  y)  which  satisfies  the  momentum  equation  and  some  of  the 
boundary  conditions,  in  a  case  where  the  velocity  profile  exponentially  decays  through 
the  boundary  layer.  As  Rosenhead  states,  “It  is  hoped  that  this  form  will  approximate 
to  the  exact  profile,  which  satisfies  all  of  the  conditions”  [Ref.  29].  The  complimen¬ 
tary  error  function  which  is  introduced  by  Timman  produces  the  desired  asymptotic 
behavior;  namely,  the  velocity  exponentially  decays  as  the  distance  from  the  surface 
increases.  We  will  attempt  to  use  a  modification  of  Timman’s  method  in  the  case  of 
the  velocity  flow  along  the  free  surface. 

2.  Timman’s  Method 

The  velocity  profile  form  assumed  is 

§  =  /M.  ”  =  ski  <V1> 

where  <5(X)  is  the  effective  total  thickness  of  the  boundary  layer  and  U(X)  is  the 
velocity  of  the  external  flow.  The  function  /  may  also  depend  on  X  through  some  of 
the  coefficients  which  are  chosen  to  be  part  of  the  profile  in  order  to  satisfy  surface 
boundary  conditions.  Timman’s  profile  is  of  the  form 

jj  =  f{r])  =  \-j  e_??2 (a  +  crj2  -j - )  dy  —  e_7?2 (6  +  dy2  -| - ),  (V.2) 

where  a,  b,  c,  and  d,  etc.,  are  coefficients  which  will  turn  out  to  be  functions  of  X. 
As  y  — >  oo,  it  follows  that  /  — ►  1  to  match  the  external  flow.  Therefore,  only  the 
boundary  conditions  at  the  no-slip  surface  rj  =  0  are  considered.  Only  the  same 
number  of  boundary  conditions  as  unknown  coefficients  are  used. 

In  Timman’s  profiles,  the  velocity  difference  f(rj)  —  1  =  ^  —  1  exponentially 
decays  through  the  boundary  layer  and  satisfies  u  =  0  at  the  surface.  In  the  core 
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region,  77  — >  00,  f(rj)  — >  1,  and  u  — >  U,  or  the  velocity  in  the  boundary  layer  matches 
the  core  velocity. 

In  our  problem,  the  surface  is  a  free  surface,  and  the  no-slip  condition  of  Tim- 
man’s  method  does  not  apply,  since  0  at  the  surface.  The  velocity  in  our  modifi¬ 
cation  must  be  adjusted  accordingly  to  match  the  free  surface  boundary  conditions. 
The  velocity  u  stills  decays  exponentially  through  the  boundary  layer,  but  where  Tim- 
man’s  approach  assumes  ^  =  0  at  the  surface  and  ^  =  1  at  the  core/boundary-layer 
edge  (where  the  boundary  layer  meets  the  core  flow),  our  problem  must  be  modeled 
to  have  u  as  a  maximum  at  the  surface  and  u—U  —*  0  at  the  core/boundary-layer 
edge,  since  from  Canright’s  numerical  data  [Ref.  1]  and  the  scaling  analysis,  the  core 
velocity  is  negligible. 

Consider  again  the  momentum  equation  in  the  free  surface  boundary  layer: 

UUX  +  Vuy  —  Uyy  (V-3) 

with  boundary  conditions  of  V  =  0  and  uy  =  Tx  when  y  —  0.  Combine  u  times  the 
continuity  equation  with  this  momentum  equation: 

U  ( UX  +  Vy  =  0)  ©  (UUX  +  Vuy  =  Uyy) 

which  gives 

UUX  +  UVy  +  UUX  +  VUy  =  Uyy 

This  can  be  rewritten  as 

(UU)X  +  (vV)y  =  (Uy)y.  ( V-4) 

Integrate  this  expression  with  respect  to  y  and  simplify: 

rOO  pOO  pOO 

/  (uu)xdy+  /  (uV)ydy  =  /  (uy)ydy  (V.5) 

Jo  Jo  Jo 

poc  1 00  1 00 

2  Jo  uuxdy  +  uV\ 0  =  wy|0  =  —  (uy)y=o  (V-6) 

Since  u  =  V  =  0asy— >00  and  V  =  0  at  y  =  0,  the  expression  for  uy  along  the 
surface  becomes 

(uy)y=o  =  ~2[  uux  dy.  (V.T) 

JO 
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To  begin  the  approximation,  use  a  value  for  Tx  =  uy,  the  free  surface  boundary 
condition.  Along  the  surface,  Tx  decreases  sharply  in  magnitude  as  the  distance  from 
the  wall  increases,  decaying  along  the  length  of  the  boundary  layer.  This  can  be  seen 
in  Figure  4,  which  shows  the  gradient  of  the  surface  temperature  in  the  r-direction, 
using  Canright’s  numerical  temperature  data.  It  also  shows  the  temperature  profile 
(dotted  line)  along  the  surface. 

Can  the  revised  momentum  equation  be  put  into  a  form  that  has  a  displace¬ 
ment  and  momentum  thickness?  Expand  the  upper  limit  of  integration  in  the  mo¬ 
mentum  equation  to  infinity,  to  make  it  consistent  with  the  range  of  u(X,  y).  Define 
82  as  the  momentum  thickness  of  the  boundary  layer.  Let 

S2  =  u2  dy  (V.8) 

Since  77  =  |,  where  8(X)  is  the  surface  boundary-layer  thickness,  y  =  8rj,  and  dy  = 
8  dr). 

If  }(rj)  is  defined  by  f(ij)  =  u,  then 

82  =  8  f  f2(rj)  dr).  (V.9) 

J  0 

Thus,  the  skin  friction,  or  shear  stress,  can  be  written  in  terms  of  the  momentum 
thickness,  yielding  a  first-order  differential  equation  r0  =  —  ^  (<52),  or 

JL(82)  +  t0  =  0  (V.10) 

3.  Modification  to  Timman’s  Method 

Consider  a  family  of  velocity  profiles  where  all  of  the  coefficients  are  functions 
of  X,  as  shown  in  (V.2).  If  no  terms  of  higher  order  than  the  c(X)  and  d(X)  terms 
are  included,  then  the  modification  to  Timman’s  profile  yields 

/*oo 

u(X,  y)  =  f(r))  =  J  e~t2  (a(X)  +  c{X)t2)  dt  -  e^2  (1 b(X )  +  d(X)r)2)  ,  (V.ll) 
where  77  =  s^x)-  When  y  =  0,  the  surface  velocity  becomes 

u(X,  0)  =  ^a(X)  +  ^c(X)  -  b(X).  (V.12) 
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Tx 


Figure  4.  Temperature  Gradient  along  the  Free  Surface 


Due  to  the  complications  of  this  model,  we  simplify  it  even  further  (see  Ap¬ 
pendix  A  for  our  attempt  at  deriving  a  full  Timman-type  model):  let  b(X )  =  c(X)  = 
d(X)  =  0.  Then 

fOG  „ 

U  =  f(rf)  =  /  a(X)e~t  dt,  (V.13) 

Jr) 

or 

f(v)  =  [1  -  erf  (77)]  =  erfc(77) ,  (V .  14) 

where  erfc(77)  is  the  complementary  error  function  of  77.  A  typical  characteristic  of 
the  complementary  error  function  of  77  is  that  as  77  increases,  erfc(77)  quickly  decays 
to  0. 

Solving  the  ordinary  differential  equation  (V.10),  (62)  +  To  =  0,  for  the 

momentum  thickness  yields 

62  =  -  [X rodX'  +  d  (V.15) 

Jo 
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However,  the  shear  stress  tq  =  uy  =  Tx  from  the  given  surface  boundary  condition, 
so 


62 


[XTX'dX'  +  cu 
Jo 


or  62  =  —  T  +  ki.  The  boundary  condition  that  must  be  satisfied  is  that  as  X  — »  00, 
82  (X)  — »  0,  u  — >  0,  T  — 0,  and  tq  — »  0.  This  implies  that  k\  =  0.  Thus, 


62  =  -T.  (V.16) 

(Since  the  temperature  is  less  than  zero,  82  is  a  positive  quantity).  The  momentum 
thickness  has  a  value  equal  to  the  magnitude  of  the  temperature.  This  implies  two 
things.  First,  the  choice  for  the  velocity  profile  does  not  explicitly  affect  the  solution 
of  the  momentum  thickness  integral.  Second,  since  the  temperature  is  of  order  1, 
the  momentum  thickness  is  of  order  1.  The  question  that  remains  is:  what  is  the 
boundary-layer  thickness  81 

Since  f(rf)  =  \(^{X)  y/7rerfc(^),  substitute  this  expression  into  the  equation 
for  the  momentum  thickness  (V.9)  and  obtain 


82  =  -T 


(V.17) 


- 


[erfc(77)]2  dr] 


(V.18) 


Now  determine  8(X)  in  terms  of  the  temperature.  First,  the  coefficient  a(x) 
needs  to  be  determined.  Recall  that  f(rj)  =  u,  so  that 


u 


y  ~ 


df 

dy 


1 

8(X) 


where  the  prime  denotes  differentiation  with  respect  to  77.  Thus,  f'(f])  =  8(X)  uy  = 
8(X)Tx  along  the  surface. 

Differentiating  the  velocity  profile  yields 


/'(*?)  =  ~a(X)  e~< 


(V.19) 
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Along  the  surface,  y  =  r\  =  0,  so  that  /'( 0)  =  —a(X),  which  implies  a(X)  = 
-6(X)  uy  =  -6(X)  Tx.  Thus 

rp2  7 r 

S2  =  ^r-  63{X)  /  [erfc(7/)]2  dy  =  — T  (V.20) 

4  Jo 

This  allows  the  determination  of  the  boundary-layer  thickness  S(X)  to  be 

rp 

HX)=  -T2-  ■  (V.21) 

.  1X£. 

where 

£=i  L  ferfc^^2  dr] = (:  _  ^r)  •  ^V-22^ 

Numerically,  e  «  0.2596  (see  Appendix  B  for  derivation).  So,  6(X)  is  positive  and  of 
order  1,  which  is  not  unreasonable. 

Assembling  the  components  of  the  profile,  the  velocity  can  be  modeled  as 

u  =  -~£~v 3  (-T)1/3  (Tx)1/3  erfc(7 ?).  (V.23) 

Examining  each  term,  (— T)1//3  decays  to  0  as  the  distance  from  the  wall,  X, 
gets  large.  (Tx)1/3  also  decays  to  0  as  X  increases.  In  addition,  erfc (rj)  decays 
rapidly  to  0  as  the  distance  from  the  surface,  y,  increases,  which  shows  that  u  goes 
to  0  at  the  edge  of  the  surface  boundary  layer.  As  seen  in  Figure  5,  there  is  excellent 
agreement  between  the  predicted  velocity  profile  along  the  surface  (discrete  points) 
and  the  velocity  profile  generated  from  Canright’s  numerical  data  (solid  curve).  In 
order  to  evaluate  the  velocity  u,  Canright’s  surface  temperature  data  was  used.  Also, 
taking  a  vertical  slice  of  u  through  the  boundary  layer  also  gives  reasonably  good 
agreement  between  the  predicted  profile  and  the  numerical  profile.  Typical  profiles 
for  a  y-direction  slice  of  the  horizontal  component  of  the  velocity,  u,  are  shown  in 
Figure  6.  As  can  be  seen,  the  value  of  u  increases  exponentially  as  the  depth  below 
the  surface  decreases.  In  this  figure,  u  increases  by  more  than  an  order  of  magnitude 
as  it  passes  through  the  boundary  layer.  The  predicted  velocity  is  given  at  distances 
of  X  =  0.0299  and  X  =  0.0494  (dotted  fines)  from  the  wall,  and  are  compared  to 
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Surface  Velocity  Comparison 


Figure  5.  Comparison  of  Surface  Velocity  Profiles  using  Timman’s  Method  Prediction 
and  Canright’s  Numerical  Data 

velocity  profiles  taken  at  the  same  ^-distances  using  Canright’s  data  (solid  lines).  The 
profiles  using  Timman’s  model  decay  much  faster  than  the  profiles  of  the  numerical 
data  as  the  distance  from  the  surface  increases.  This  might  be  due  to  the  lack  of 
a  second  parameter  in  the  model.  Also,  note  that  the  surface  boundary  condition 
(uy  ~  Tx)  matches  for  each  corresponding  plot. 

B.  VISCOUS  WALL  BOUNDARY  LAYER 
1.  Background 

The  behavior  of  the  velocity  along  the  vertical  rigid  wall  is  similar  to  that 
outlined  in  Timman’s  paper.  The  no-slip  condition  exists  at  the  wall,  so  the  vertical 
component  of  velocity  is  zero  at  the  wall,  increases  to  a  peak  position  (where  the 
vorticity  is  zero),  and  then  decreases  to  match  the  value  of  the  core  velocity  at  the 
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Figure  6.  Typical  Profiles  of  Velocity  (u)  through  the  Surface  Boundary  Layer 

edge  of  the  viscous  boundary  layer.  Two  approaches  are  studied.  The  first  is  a  method 
by  Glauert,  which  considers  the  velocity  profile  through  the  wall  boundary  layer  as 
similar  to  that  near  a  plane  wall  jet.  The  Glauert  approximation  matches  numerical 
work  done  by  Canright  [Ref.  1]  in  1994  quite  well.  A  discussion  of  the  vorticity  flux 
is  also  included  as  an  appendix,  as  the  vorticity  u  can  be  represented  by  the  gradient 
of  the  vertical  component  of  the  velocity,  vx.  An  interesting  conclusion  is  reached. 
The  second  method  is  a  Timman-type  approximation,  but  since  four  parameters 
are  needed  for  accuracy,  this  method  involves  solving  coupled  nonlinear  differential 
equations.  Since  the  Timman  method  worked  with  the  surface  boundary  layer,  its 
use  with  the  wall  boundary  layer  was  investigated  and  is  given  as  an  appendix.  The 
Glauert  approach  is  more  desirable  due  to  its  simplicity. 
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2.  Glauert’s  Method 

Suppose  the  fluid  flow  from  the  free  surface  into  the  cold  corner  turns  down  the 
wall  and  acts  in  the  manner  of  a  plane  wall  jet.  There  is  a  similarity  solution  for  the 
velocity  profile  through  the  wall  boundary  layer.  Following  Glauert  [Ref.  36],  start 
with  the  dimensionless  momentum  equation,  assuming  the  pressure  is  everywhere 
uniform, 

U  vx  +  v  vY  =  vxx.  (V.24) 

Combine  the  continuity  equation  with  a  definition  of  a  stream  function  ip,  defined 
by  U  =  §  and  v  —  —  ff.  Consider  the  possibility  that  there  shall  be  a  similarity 
solution  of  these  equations,  with  the  velocity  proportional  to  Ya  and  the  jet  thickness 
proportional  to  Yb.  The  two  sides  of  the  momentum  equation  (V.24)  will  vary  with 
Y  in  the  same  manner  if  a  +  26  =  1.  For  a  wall  jet,  fluid  momentum  is  not  conserved, 
due  to  shear  stress  from  the  wall. 

Write  the  stream  function  as 


=  -Y'-bf{r,) 


where 


Then 


rj  =  (1  -  b)xY~b. 


v  =  (1  - 

and  (V.24)  becomes 

/"'  +  ff"  +  af2  =  0  (V.25) 

where  a  =  (26  —  1)/(1  —  6).  The  boundary  conditions  require  that  /( 0)  =  /'( 0)  =  0 
and  /'( oo)  =  0.  Integrating  (V.25)  between  the  limits  r)  and  oo  yields 

roc 

f  +  ff  -  («  -  1)  /  fndr)  =  0.  (V.26) 

Jr} 

Let  g(rj)  =  fjj°  f2dr),  such  that 

f  +  ff-(a-l)g  =  0.  (V.27) 
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Notice  that  g  >  0.  Multiply  both  sides  of  the  above  expression  by  /',  yielding 


/7"  +  //'2- (a -!)»/'  =  0.  (V.28) 


Integrate  again,  using  integration  by  parts  on  the  last  two  terms  and  thus 

i/'2  -/»  +  (<*-  2)  jT  S' 9  *)  =  0.  (V.29) 

If  /  =  0  when  /'  =  0,  then  a  =  2,  since  /“  f'gdr)  >  0.  This  gives  a  =  — |  and  6  = 

Regarding  what  Glauert  terms  the  external  momentum  flux,  start  again  with 
the  momentum  equation  (V.24).  Integrate  this  with  respect  to  x  between  the  limits 
x  and  oo,  using  the  conditions  that  v  — ►  0  as  x  — *■  oo. 


U  vxdx  + 


vvy  dx 


vxx  dx 


(V.30) 


By  continuity,  Ux  =  —vy,  and  from  integration  by  parts,  the  first  term  above  becomes 

roo  roc 

/  U  vxdx  —  —Uv  +  /  vvydx  (V.31) 

Jx  Jx 


Substituting, 


d  f°°  o 

—  y  v2dx-  Uv  +  vx  =  0  (V.32) 

Multiply  this  equation  by  v  and  integrate  with  respect  to  x  between  0  and  oo. 

f)  roc  /  roo  \  roc  /  roc  \  roc  1  |OQ 

W I  V(L  v2dX)  dX~J  Vy  ( L  y2dX)  dX  ~  Jo  y2U  dX  +  2V1o  =  °  (V-33) 

The  last  term  on  the  left  is  zero,  and  from  continuity, 

roo  /  roo  \  roo  ^  IOO  roo  „ 

—  J  Vy  IJ  v2dx'\dx  =  U  j  v2dx\^-\-  j  v2U dx 


Thus, 


0 

dY 


dx  =  0, 


(V.34) 
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since  U{ 0)  =  w(0)  =  0. 

Let 

F  =  J  v  (^J  v 2  dx'^j  dx  (V.35) 

Then  =  0,  and  F  =  constant,  independent  of  the  distance  Y  from  the  free  surface. 
However,  realistically,  there  may  be  some  interference  from  the  velocity  leaving  the 
surface  and  coming  from  the  corner.  Glauert  states  that  physically,  the  “flux  of  the 
exterior  momentum  flux  is  constant,  but  this  is  hardly  a  familiar  concept.”  [Ref.  36]. 

It  is  worth  pausing  here  to  find  a  rough  estimate  of  the  magnitude  of  the  flux  of 
the  exterior  momentum  flux,  F.  With  a  knowledge  of  the  conditions  in  the  impinging 
free  jet,  based  on  the  entry  conditions  from  the  corner,  start  with 


F  =  Io°V  0f°  ^  dX' )  dX  ~V*  fo°V  (/°°  w  dx>)  dx’  (V.36) 

where  v*  is  a  typical  jet  velocity.  Write 


UM-HM 

(V.37) 

So, 

1  r°°  d  /  \ 2 

-aid vix) dx 

(V.38) 

However, 

which  concludes  that 

F*\v’{Cvdxf 

(V.39) 

This  means  that 

F  =  ^(typical  velocity)  x  (mass  flux)2 
z 

(V.40) 

in  the  plane  jet  case.  Neither  the  mass  flux  nor  the  magnitude  of  the  fluid  velocity 
will  vary  much  when  a  plane  jet  is  deflected  off  of  the  vertical  wall,  so  (V.40)  is  a 
good  estimate  of  the  constant  flux  of  exterior  momentum  flux,  F.  A  more  accurate 
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calculation  of  F  will  be  given  later.  Keep  in  mind  that  this  approximation  for  F  has 
a  singularity  at  the  origin,  where  Y  =  0. 

The  similarity  form  of  the  velocity  distribution  depends  only  on  F  and  is  given 
Glauert’s  paper  [Ref.  36]  as 

/"'  +  //"  +  2/'2  =  0  (V.41) 


with  boundary  conditions  /( 0)  =  /'( 0)  =  0  and  /( oo)  =  0.  (As  an  aside,  Batchelor 
refers  to  a  steady  narrow  two-dimensional  jet  of  fluid  adjoining  a  plane  rigid  wall  and 
offers  4 /'"  +  //"  +  2 /'2  =  0  as  the  differential  equation  [Ref.  37].) 

Start  with  Glauert’s  equation,  (V.41).  Multiply  the  equation  by  /.  Integrate 
the  first  term  by  parts: 

/  s r" = //"  -  Jf'f" = //"  -  If2 

The  second  and  third  terms  in  the  velocity  distribution  form  a  perfect  differential, 

/  (ff + vs'2)  =  s2f 


Thus, 

//"  ~ +  ff  =  const  =  0,  (V.42) 

since  /'(oo )  =  0. 

Multiply  now  by  /“V2: 

ri/2r  -  \rz,2f2 + /1/2r = o.  (v.43) 

Integrate  again  to  obtain  (this  time,  the  first  two  terms  form  a  perfect  differential) 

f~1/2f  +  |/3/2  =  constant  (V.44) 

o 

If  fo(ri)  is  a  solution  of  the  original  third-order  differential  equation,  then  so 
is  fi(rj)  =  Afo(Arj)  for  any  constant  A,  and  /i  satisfies  the  boundary  conditions  if  /o 
does.  According  to  Glauert  [Ref.  36],  without  loss  of  generality,  select  the  solution 
with  /(oo)  =  1  and  take  the  above  constant  to  have  the  value  of 
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Now,  write  /  =  g2.  Then  f  =  2 gg'  and  the  differential  equation  in  /  becomes 

=  (V.45) 

which  gives  as  a  solution 

(V.46) 


Figure  7.  Velocity  Profile  from  Glauert’s  Method 

A  plot  of  the  variation  of  Glauert’s  velocity  profile  (/'(tj)  vs  77 )  is  given  in 
Figure  7.  This  profile  is  similar  to  the  profile  determined  using  Timman’s  approach 
and  to  the  profile  obtained  using  the  numerical  data  from  Canright  [Ref.  lj. 

The  velocity  and  position  equations  may  be  'written  as 


(§f/» 

(V.47) 

( 5F  y/4* 

V32W 

(V.48) 
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and 


Recalling  that  77  =  the  above  equation  can  be  used  to  give  an  expression  for  the 
boundary-layer  thickness,  A: 


A  = 


5  F  \~1/4 
32  Y3J 


(V.49) 


There  is  a  drawback  to  Glauert’s  method.  The  momentum  flux  has  a  singular¬ 
ity  when  Y  =  0,  which  will  not  allow  the  velocity  into  the  corner  (along  the  surface 
at  X  =  0)  to  match  the  velocity  out  of  the  comer  (along  the  wall  at  Y  =  0).  This 
means  that  the  surface  velocity  is  based  on  one  ratio  of  the  temperature  gradient 
(u  oc  Tx3),  while  Glauert’s  suggested  matching  (V.40)  bases  the  wall  velocity  on 
the  inverse  ratio  (v  oc  7^1/3).  Thus,  along  the  surface,  as  the  temperature  profile  is 
compressed,  the  velocity  into  the  corner  increases.  However,  along  the  wall,  as  the 
isotherms  are  compressed,  the  velocity  decreases,  which  is  not  reasonable. 

As  an  alternative  matching  condition,  to  eliminate  the  singularity  of  the  simi¬ 
larity  solution  at  the  origin,  calculate  F  when  Y  is  a  distance  equal  to  the  boundary- 
layer  thickness,  A.  This  will  give  a  strength  value  for  F  which  is  relatively  equal  to 
the  value  at  the  comer,  but  the  singularity  is  now  avoided.  Let 


Solving  for  Y  yields 


5  F 


(V.51) 


Now,  matching  the  magnitude  of  the  surface  velocity  expression  into  the  corner  to 
the  maximum  wall  velocity  expression  out  of  the  corner,  we  find 


f  (ff=2-5/3( 


/5  F\1/2 
\2Y  ) 

(V.52) 

1/3 

1 

(V.53) 

which  reduces  to 


This  matching  condition  for  F  embodies  the  positive  feedback  mechanism  of  the  cold 
corner. 
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C.  CORE  FLOW  REGION 

I  NEVER  COME  ACROSS  ONE  OF  LAPLACE’S  “THUS  IT  PLAINLY  APPEARS'  WITHOUT  FEEL¬ 
ING  SURE  THAT  I  HAVE  HOURS  OF  HARD  WORK  BEFORE  ME  TO  FILL  UP  THE  CHASM  AND 
FIND  OUT  AND  SHOW  HOW  IT  PLAINLY  APPEARS.  — N.  Bowditch 


1.  The  Green’s  Function 

Recall  that  the  core  velocity  U  is  treated  as  incompressible  irrotational  flow. 
As  explained  earlier,  a  stream  function  ip  is  introduced  such  that  V2ip  =  0.  Define 
a  Green’s  function  G(x,x o)  for  Laplace’s  equation.  Thus,  V2G(x,  xo)  =  S(x  —  xq), 
where  x0  is  an  arbitrary  point  in  the  domain.  The  solution  can  then  be  represented 
in  integral  form.  First,  however,  the  boundary  conditions  in  the  core  flow  region 
must  be  prescribed  in  terms  of  ip.  The  core  flow  must  match  the  entrainment  into 
the  boundary  layers,  that  is,  the  normal  velocities.  Along  the  wall  boundary  layer, 
U(x  — ►  oo,  Y)  =  ipy{X  — »  0,  Y).  Along  the  surface  boundary  layer,  V(X,y  — ►  oo)  = 
—ipx(X,  Y  —>  0).  These  conditions  can  be  simplified  to  Dirichlet  conditions,  matching 
the  stream  function  ip  at  the  outer  edges  of  the  boundary  layers. 

The  Green’s  function  is  determined  using  the  method  of  images.  Start  with 
the  solution  in  a  quarter  plane  and  use  symmetry  in  the  other  quadrants  to  match 
the  boundary  conditions.  There  axe  four  quadrants  with  logarithmic  singularities. 
The  defining  problem  for  the  Green’s  function,  V2G(x,  x0)  =  8(x  —  x0),  satisfies  the 
corresponding  homogeneous  boundary  conditions,  G(x,  0;  xo,  yo)  =  G( 0,  y;  x0,  y0)  =  0. 
Here,  the  notation  that  G(x,x o)  =  G(x,y,x0,yo)  is  used.  The  quarter-infinite  space 
(x  >  0,  y  >  0)  has  no  sources  except  at  a  concentrated  source  at  x  =  xo. 

Using  the  method  of  images  in  the  core  flow  region,  the  Green’s  function  for 
Laplace’s  Equation  on  a  quarter-plane  (X  >  0,  Y  >  0)  is 

G(x,x0)  =  ^{in[(x-xo)2  +  (y-y0)2]-in[(jf-xo)2  +  (y+n)2] 

-In  [(X  +  X0)2  +  (Y-  V0)2]  +  In  [(X  +  X0f  +  (K  +  F0)2] } 
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-  1  ■  f  K*  -  X°)2  +  (Y  -  n)2)  [(X  +  Xcf  +  <Y  +  Vo)2]  1  , 

^  U(x-x0r+(Y+Y0)m(x+x^+(Y-Y0nS  ' 

To  check  the  boundary  conditions,  start  with  the  Green’s  function  and  set  X 
or  Y,  respectively,  equal  to  zero.  Thus  it  plainly  appears  that  G  =  C?  =  0. 
Since  the  boundary  condition  for  the  Green’s  function  is  the  Dirichlet  kind  at  Y  —  0, 
and  at  X  =  0,  then  a  positive  image  source  must  be  used  for  Y  <0  and  X  <  0,  and 
negative  sources  must  be  used  for  X  >  0,  Y  <0  and  for  X  <  0,  Y  >  0.  In  this  way, 
equal  but  opposite  sources  are  located  at  the  symmetry  locations  of  Xo  in  the  four 
quadrants  of  the  plane,  and  there  is  no  flow  along  Y  —  0  and  X  =  0,  as  desired. 

To  solve  Laplace’s  equation,  multiply  the  Green’s  function  identity  by  ip,  giving 

[v2G(X,Xo)  =  6(X-Xo)]. 

Multiply  Laplace’s  equation  by  G,  yielding 

G  [v2ip  =  0] . 

Subtracting  and  integrating  over  the  domain  A  gives 

J  V2G  -  G  VV)  dA  =  i>  (X0)  (V.55) 

Using  Green’s  second  identity,  Laplace’s  equation  in  two  dimensions  becomes: 

J  Ja  (i>V2G  -  G  V2V>)  dA  =  jf  (GVxp  -ijfVG)  •nds  =  'ip  (X0)  ,  (V.56) 

where  n  is  the  normal  to  the  boundary  (inward  normal)  and  s  is  the  arc  length  along 
the  boundary.  Consider  a  large  quarter-circle  (in  the  limit  as  the  radius  tends  to 
infinity).  The  inward  normal  for  the  top  (boundary  layer  along  the  free  surface)  is  j 
and  the  inward  normal  for  the  side  (boundary  layer  along  the  wall)  is  i.  It  will  be 
shown  that  the  contribution  at  oo  tends  to  vanish  if  ?/>  — >  0  sufficiently  fast. 
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(V.57) 


+  lim  l{GXiP-ipXG)-hds 

|X|— oo  Ja 

The  last  term  on  the  right  is  the  contribution  as  the  radius  of  the  quarter  circle  (or 
|X|)  tends  to  infinity.  It  will  vanish  if  -0  — s-  0  sufficiently  fast,  or  if  ^  ~  To  show 
this,  integrate  in  polar  coordinates,  for  which  ds  =  rdd  (r  =  |X|): 


hm  f  (GV^  -xpXG)  rdO  =  0. 

Substituting  G  and  evaluating  /  dd  =  2w  yields  the  necessary  condition 

/  dib  \ 
hm  r  lnr  — - yj\  =0. 


(V.58) 


dr 


(V.59) 


If  ^  0  sufficiently  fast,  or  if  ^  ~  then  the  contributions  of  this  “boundary” 

\X\ 

term  wifi  vanish. 

Hence, 

dip 


dG 

dX. 


dY  + 


x=o 


r{ew-*w)  L"  <v-6o) 


Since  G 


x=o 


=  G\=  0, 


y= o 


*w-f(-+ff)  (-*§?) 


lx=0 

and  V  =  -U 

G  —  0atX  =  Y  =  0,  and  -  J| 


dX 


Y= 0 


(V.61) 


Recall  that  U  =  U 

dY  x=o 


dX 


y=o 


.  Since  the  boundary  conditions  give 


=  0,  we  only  need  to  find  dGXx°)  which  is 
x=o  ’  J  dY  y=o’ 


(V.62) 


dG(X,X0) 

_  1 

n  Y0 

dY 

y=o—  t r 

[(X  +  X0)2  +  r02  (X  -  x0)2  +  Y02\ 

Thus,  for  the  core  flow  due  to  the  surface  boundary  layer,  insert  the  mass  flux  into 
the  corner  for  yielding 

*  =  -/% W«0)3G(1’J?o)l 

Jo 


dY 


y= o 


dX, 


(V.63) 


As  stated  previously,  one  advantage  of  the  Green’s  function  is  that  the  flow 
can  still  be  represented  over  the  entire  quarter  plane,  so  that  there  is  no  artificial 
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recirculation  due  to  imposed  artificial  boundaries.  Thus,  only  the  boundary  condi¬ 
tions  at  the  boundary  layer/core  region  interfaces  are  required.  The  Green’s  function 
approach  does  not  require  flow  boundary  conditions  at  the  computational  domain’s 
boundaries. 


2.  Mass  Flux  Calculations 

In  order  to  obtain  a  solution  to  Laplace’s  Equation  in  the  core  flow  region,  the 
boundary  conditions  from  the  surface  and  wall  boundary  layers  must  be  addressed. 
Recall  that  along  the  surface  boundary  layer,  the  velocity  is  modeled  as 

u  =  ^  e~1'3  (-T)1/3  ( Tx erfcfr).  (V.64) 


As  explained  previously,  the  horizontal  velocity  component  u  decays  to  zero  as  X  and 
rj  increase. 

To  determine  the  mass  flux  along  the  surface  layer  into  the  cold  corner,  sum  the 
velocity  component  as  (the  correct  scalings  for  the  velocity  and  the  boundary-layer 
thickness  are  inserted  as  appropriate) 


^surface 


(_ T)i/3  (Tx)^edc(V) 


dy 


=  -JQ  l£  1/3(— r)1/3(r*)1/3erfc(r7)]  6  dy 


=  [e-1^  (-10*^ 


=  s~2/H-T)2/3  (Tx)-1'3  jTerfcO,)*, 


(V.65) 


e“2/3  (— T)2/3  (Tx)~1/3 

A 
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When  X  =  0,  the  stream  function  (mass  flux)  into  the  corner  can  be  calculated  as 
\  £~2^  (—To)2/3  (TXo)~1/3.  This  value  will  be  compared  to  the  mass  flux  out  of  the 
corner  shortly. 

To  determine  the  mass  flux  down  the  wall,  sum  the  vertical  velocity  compo¬ 
nent,  v.  Recall  that 


where 

x  /  5  F  \1/4 
^  =  A  =  \32yv  *• 

Thus,  the  stream  function  (mass  flux  with  appropriate  scalings  for  the  velocity  and 
boundary-layer  thickness)  is 


=  -(5F)>'4(2)3',4r1/4/{>?)|”  (V.66) 

We  have  shown  that  /( 0)  =  0  and  lim^^oo  /(r?)  =  1,  so  that  the  mass  flux  down  the 
wall  is 

Vwn  =  -(40Fy)1/4  (V.67) 


D.  POTENTIAL  FLOW  DUE  TO  A  WALL  JET 

The  next  test  for  our  problem  is  to  check  the  potential  flow  due  to  a  wall 
jet  in  the  quarter-plane.  A  jet  is  a  flow  in  which  the  width,  or  cross-stream  scale,  is 
considerably  smaller  than  the  downstream  scale.  For  a  wall  jet,  this  type  of  flow  occurs 
along  a  solid  boundary.  In  Glauert’s  discussion  of  the  wall  jet,  he  postulated  that  the 
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entire  flow  of  the  wall  jet  cannot  conform  to  one  similarity  solution  over  the  entire 
domain  [Ref.  36].  Glauert  divided  the  flow  region  into  inner  and  outer  portions, 
treating  the  two  portions  separately.  The  dividing  line  for  the  two  portions  was 
the  point  of  maximum  velocity.  Additionally,  as  we  have  demonstrated  previously, 
Glauert  also  showed  that  the  “flux  of  exterior  momentum  flux”  across  the  viscous 
boundary  layer  is  constant  [Ref.  36] .  The  plane  wall  jet  is  created  when  a  jet  flow 
is  discharged  in  a  tangential  direction  along  the  plane  surface;  in  this  case,  the  wall. 
Plotkin  [Ref.  40]  studied  a  higher-order  correction  to  Glauert ’s  solution,  providing  a 
second-order  solution  including  displacement.  The  first-order  inner  viscous  solution 
induces  a  correction  in  the  outer  inviscid  core  flow  solution.  Plotkin  matched  the  inner 
and  outer  solutions  and  developed  a  solution  which  satifies  the  boundary  conditions 
in  a  semi-infinite  plane,  where  ip(x,  0)  =  4x1//4  for  x  >  0,  and  ^  =  0  for  x  <  0.  His 
solution  satisfies  Laplace’s  equation  and  is 


ip(x,  y)  =  4  ^R.(x  +  iy )1^  —  S(a:  +  iy)1/4]  (V.68) 


where  K  and  S  denote  the  real  and  imaginary  parts,  respectively,  for  0  <  arg(x-Hy)  < 
27 r. 


In  the  present  problem,  Plotldn’s  solution  must  be  modified  to  match  the 
boundary  conditions  in  the  quarter-plane,  namely,  that  ^(0,  Y)  =  — (40FY)1/4  and 
when  Y  =  0  and  X  >  0,  V’  =  0.  The  angle  is  |,  and  after  taking  the  fourth-root,  the 
argument  for  the  complex  variable  becomes  |.  Taking  into  account  that  the  jet  flows 
down  the  wall  (in  the  y-direction),  the  potential  solution  becomes 


$(X,Y)  =  —(40  F)1/4 


K{Y  +  iX)1/A  - 


V2  +  vga(y+,X),/4 

y  2  —  y/2 


(V.69) 


which  simplifies  to 


1>(X,  Y )  =  -(40F)1/4  [K{Y  +  iX)1/4  -(1  +  V2)  3(Y  +  iX)l/A]  .  (V.70) 


This  also  satisfies  Laplace’s  equation  in  the  quarter-plane.  This  expression  for  ip(X,  Y) 
is  used  to  determine  the  core  velocity  field  due  to  the  wall  boundary  layer. 
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VI.  NUMERICAL  RESULTS 


REFINEMENT  OF  THE  METHOD  ALSO  REQUIRES  A  MORE  REFINED  USER.  —  Federico  Paris 
and  Jose  Canas 


The  PURPOSE  OF  COMPUTATION  IS  INSIGHT,  NOT  NUMBERS.  —  Richard  W.  Hamming 


A.  SOLVING  LAPLACE’S  EQUATION 

Before  solving  Laplace’s  equation  in  the  core-flow  region,  the  numerical  code  to 
solve  Laplace’s  Equation  was  tested  on  a  simple  problem  with  a  known  solution,  the 
point  vortex.  The  two-dimensional  solenoidal  flow  in  the  core  flow  region  associated 
with  a  straight  line  vortex  may  be  described  in  terms  of  a  stream  function  as 


^  =  -Al°g(T,  (VL1) 

where  k  is  the  strength  of  circulation  and  a  =  \J{x  —  £0)2  +  (y  —  Vo)2-  In  a  flow  field 
which  is  entirely  two-dimensional,  the  appropriate  term  for  the  singularity  is  point 
vortex.  The  line  vortex,  a  line  in  space  that  has  a  given  circulation,  or  vorticity, 
becomes  a  point  vortex  when  a  single  plane  cuts  through  it.  For  simplicity,  set  k  =  1 
and  find  0,  the  harmonic  conjugate  of  ip.  Since  the  two  scalar  functions  <p(x,  y)  and 
ip(x,  y)  both  satisfy  Laplace’s  equation,  they  also  provide  alternative  specifications  of 
a  core  flow  vector  U  which  is  both  irrotational  and  solenoidal  [Ref.  37].  Thus,  given 
ip  above,  and  the  fact  that 


dip 

dy 


d(p 

dx 


(p  is  found  to  be 


K, 

2tx 


,-i 


&P  =v  _  d<P 

dx  dy  ’ 

(VI-2) 

(;) 

(VI-3) 
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To  find  the  velocity  anywhere  in  the  fluid,  find  the  distance  and  direction  to 
the  point  vortex.  The  direction  of  the  fluid  is  always  a  circle  surrounding  the  point. 
The  magnitude  of  the  velocity  is  given  by 

0  -  2 £■  (VL4) 

The  horizontal  and  vertical  components  of  velocity,  evaluated  at  x  =  0  and  y  =  0, 
respectively,  axe 

u  =  —  __I,  .1 

dx  *=o  y  ’  dy  j/=°  x 

The  Laplace’s  equation  solver  was  rim  on  this  point  vortex  problem,  treating 
the  point  vortex  at  the  origin  and  evaluating  the  <p  distribution  in  a  quarter-plane  with 
a  uniform  grid.  A  trapezoidal  scheme  for  integration  is  used.  The  output  matched  the 
behavior  of  the  known  solution,  as  seen  in  Figure  8.  The  plot  of  the  tan_1(y/x)  gives 
straight  fines  through  the  origin.  The  numerical  code  integrates  to  regions  beyond 
the  area  of  interest  (due  to  the  Green’s  function),  and  then  the  area  of  interest  -  the 
comer  region  -  can  be  analyzed  in  detail.  This  was  useful  to  test  the  code,  although 
the  stream  function  ip  will  be  used  instead  of  the  velocity  potential  <p. 


B. 


SOLVING  THE  HEAT  EQUATION 

The  full  two-dimensional  time-dependent  problem  can  be  expanded  as 


f&T  dT  3T\  d2T  d2T 

m  +  u  ax  +  v  ay )  -  ax2  +  dYr 


(VI.  6) 


The  boundary  conditions  for  the  quarter-plane  are  that  when  y  =  0,  Ty  =  0,  and 
when  x  =  0,  T  =  —  1.  Artificial  boundaries  were  introduced  to  modify  the  quarter- 
plane  problem  into  a  rectangle  problem.  Since  we  are  concerned  with  the  activity  in 
the  corner,  the  domain  was  chosen  large  enough  that  the  temperature  gradients  and 
velocity  values  decayed  to  zero  far  from  the  corner  (as  X,  Y  — >  0,  VT,  U,  V  — 0). 
Specifically,  on  the  right  boundary,  T  —  0  for  the  hot  incoming  fluid,  and  on  the 
bottom  boundary,  Ty  =  0  to  minimize  downstream  influence.  Along  the  surface  layer, 


54 


Figure  8.  Contour  Plot  of  Numerically-Solved  Phi  for  Point  Vortex 


the  heat  flux  Ty  condition  (IV.  17)  is  incorporated  into  the  scheme.  It  is  assumed  that 
the  time-dependent  solution  will  converge  to  a  steady-state  solution.  Hence,  the  initial 
conditions  can  be  chosen  arbitrarily. 

Equation  (VI.6)  was  solved  using  an  alternating  direction  implicit  (ADI)  method 
as  developed  by  Polezhaev  in  1967  [Ref.  39].  This  method  treats  each  time  step  as 
two  half  steps.  The  above  equation  can  be  substituted  into  the  scheme  as: 


TJ.  .  _  p 

Ul’3  2  AX  dx , 


STEP  1: 

,  At  ( 

1+t( 

STEP  2: 

1  At  (ir  by  c2 

+  Tr,'i2A Y~Sy 
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1h3  ““ 
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At  (y.  8y  p" 
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2AX  6x , 
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(VI-7) 

(VI-8) 
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where  the  central-difference  operators  <5  and  <52  are  defined  as: 

8x  Tij  =  Ti+ij  —  Ti-ij 
8y  T^j  =  Tij+ 1  —  Tij- 1 


(VI -9) 


&Z&  = 


rnn  _  o  T^n  I  T^n 


(VI.  10) 


*  iJ  (AX)2 

/Tnn  _  c\  rpn  ,  rpji 

=  wi  (Ay  ^  (vi-io) 

As  a  result  of  splitting  the  time  step  in  this  algorithm,  only  tridiagonal  systems 
of  linear  algebraic  equations  must  be  solved  (at  each  half-step).  This  method  is 
first-order  accurate  in  time  but  second-order  in  space,  with  a  truncation  error  of 
0[At,  ( Ax )2,  (Ay)2}  and  is  unconditionally  stable  for  the  linear  (no  convection)  case. 

The  scheme  was  coded  using  the  C  programming  language  with  a  uniform 
(Ax  =  Ay)  rectangular  grid.  The  problem’s  boundary  conditions  as  stated  earlier 
were  incorporated  into  the  scheme.  A  copy  of  this  code  is  attached  in  Appendix  E. 
To  set  the  boundary  conditions  at  each  of  the  four  sides  of  the  rectangle,  some  of 
the  coefficients  for  T^,T*j,  and  Tfj1  had  to  be  manipulated.  These  manipulations 
ensured  adherence  to  the  boundary  conditions  without  altering  the  solution  method 
inside  the  boundaries.  For  example,  second-order  difference  equations  for  the  first 
and  second  spatial  derivatives  were  modified  at  the  surface  and  at  the  bottom  of  the 


A  subroutine  tridiag.  c  is  called  in  the  execution,  which  solves  the  tridiagonal 
set  of  equations  Ax  =  b  for  the  unknown  left-hand-side  of  each  equation  in  STEP  1 
and  STEP  2  of  the  ADI  method.  In  this  case,  the  matrix  A  is  the  coefficient  matrix 
involving  velocity  and  time-derivative  terms,  b  is  the  vector  consisting  of  the  right- 
hand-side  of  each  equation  in  STEP  1  and  STEP  2  (the  known  temperature),  and  x 
is  the  vector  of  the  unknown  temperature. 

The  method  was  tested  against  known  problems  involving  mixed  boundary 
conditions  and  various  scenarios  of  velocity,  with  and  without  convective  terms,  and 
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after  iteration,  a  steady-state  was  reached,  i.e.,  the  solution  converged.  The  ADI 
method  is  unconditionally  stable  but  may  converge  slowly.  An  alternate  technique  is 
to  update  each  row  or  column  in  the  step  temperature  as  it  is  solved.  This  leads  to  a 
faster  convergence;  however,  it  is  not  clear  if  this  modified  technique  is  unconditionally 
stable. 


The  first  sample  problems  that  the  ADI  method  were  tested  against  involved 
zero  velocity  (|£  =  V2T  -  pure  diffusion,  no  convection).  The  solutions  were  linear  in 
x :  T(x,  y)  =  where  L  is  the  distance  to  the  boundary  in  the  rr-direction.  Other 
sample  problems  involved  a  non-zero  velocity  in  one  direction  +  u  •  VT  =  V2T 
-  diffusion  and  convection).  In  the  latter  case,  one  example  used  involved  a  unit 
square  domain,  with  Dirichlet  conditions  (T  =  0)  on  the  left  and  right  boundaries 
and  Neumann  conditions  on  the  bottom  ( Ty  =  0)  and  surface: 


f(x )  =  e1/2  sin(7rx)  sinh 


/  \/l  +  47 r2' 


V 


The  velocity  field  was  (U,  V)  =  (0,  —1)  in  this  example.  The  solution  to  the  problem 
is 

(Vl+47r2-l)(y-l)  (Vl+4*2+lM«-l)  ~ 

e  2  e  2 

Vl  +  47T2  —  1  \/l  +47T2  +  1 


The  ADI  algorithm  converged  to  four-decimal  accuracy  in  49  iterations,  yielding  the 
solution  shown  in  Figure  9.  The  grid  had  a  uniform  spacing  of  0.05.  The  time  step 
per  iteration  was  0.002. 


C.  THE  ADI  METHOD 

The  ADI  method  was  used  to  solve  the  two-dimensional  unsteady  heat  equa¬ 
tion,  beginning  with  an  initial  temperature  field  which  was  exponentially  decaying 
on  the  surface  and  zero  elsewhere.  This  temperature  field  was  then  iterated  in  two 
loops:  the  inner  loop  had  a  fixed  number  of  iterations  which  was  simply  the  APT 
algorithm  marching  in  time;  the  outer  loop  updated  the  velocity  profile  with  another 
fixed  number  of  iterations.  Since  the  velocity  field  is  based  upon  the  temperature 
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After  50  Iterations  -  ADI  method 


Figure  9.  Contour  Plot  of  Numerically-Solved  Temperature  for  Test  Case 

and  its  derivatives,  the  interim  temperature  solution  (after  each  set  of  inner  loop  it¬ 
erations)  was  used  to  calculate  the  velocity,  which  was  then  fed  into  the  ADI  solver 
for  another  series  of  inner  loop  iterations,  yielding  another,  more  accurate,  tempera¬ 
ture  solution.  This  process  was  repeated  a  number  of  times  equal  to  number  of  the 
outer  loop  iterations.  Once  the  difference  between  successive  temperature  solutions 
fell  below  a  certain  tolerance  (in  this  case,  0.00001),  steady-state  was  reached.  Figure 
10  shows  a  typical  steady-state  temperature  solution.  The  spatial  step  size  is  0.01  in 
both  directions. 

The  thermal  field  is  compressed  against  the  wall  by  the  velocity  into  the  corner. 
The  problem  used  20  iterations  in  each  inner  loop,  and  it  converged  in  only  12  outer 
loop  iterations.  A  closer  view  of  the  corner  is  shown  in  Figure  11,  which  also  compares 
our  approach  (on  the  left)  with  Canright’s  numerical  data  (on  the  right).  The  com¬ 
parison  is  encouraging,  particularly  since  the  model  has  no  adjustable  parameters. 
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Steady-State  Temperature  Solution 


Figure  10.  Steady  State  Temperature  Contour  Plot 


To  get  a  better  indication  of  our  method,  surface  temperatures  are  plotted,  in 
Figure  12.  There  is  excellent  agreement  of  the  temperature  gradient  from  the  corner 
into  the  core  flow  region,  well  beyond  the  boundary-layer  thickness. 


D.  UNIFORM  VELOCITY 

Following  the  boundary-layer  matching  theory  outlined  in  Bender  &  Orszag, 
[Ref.  41],  the  uniform  velocity  for  the  entire  domain  can  be  written  as 


^uniform  —  dinner  'F  U<. 


outer 


^match 


(VI.11) 


The  inner  solution  is  composed  of  the  velocity  components  in  the  two  boundary  layers, 

tanner  =  ^surface  "F  ^Avall ,  (VI.  12) 
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Figure  11.  Steady  State  Temperature  Contour  Plot  -  Comparison  with  Canright’s 
Data 


while  the  outer  solution  is  made  up  of  the  two  components  of  the  core  velocity,  due 
to  the  similarity  solution  and  the  Green’s  function: 


U, 


outer 


—  Ug reen  “1“  Us[m  Ua 


(VI.  13) 


In  the  above  expressions,  {/surface  refers  to  the  surface  boundary  layer  velocity 
component  obtained  using  the  modified  Timman  approach,  due  to  the  thermal  gradi¬ 
ent  along  the  surface,  Tx.  t/wan  refers  to  the  wall  boundary  layer  velocity  component 
obtained  using  Glauert’s  method,  fed  from  the  surface  mass  flux.  The  core  velocity 
component  Ugreen  is  a  velocity  field  obtained  using  the  Green’s  function  (from  solving 
Laplace’s  equation  in  the  core  region)  and  due  to  the  entrainment  into  the  surface 
layer,  and  the  core  component  t/sjm  is  the  velocity  field  due  to  the  similarity  solution 
of  the  potential  from  the  wall  jet  (due  to  the  Plotkin  modification).  Determine  the 
uniform  flow  in  terms  of  the  stream  function,  ip. 

Consider  first  the  surface  velocity  vector.  As  shown  earlier, 

tWfeo  =  -^£-V3(_r)>/3(rx),/3erfc(>!).  (VI.  14) 
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Steady-State  Surface  Temperature 


Figure  12.  Steady  State  Surface  Temperature  Plot 


Let  g(rj)  =  e  *2  ds  =  &  erfcfo),  where  rj  =  y/  jfr  ■  Define  G (r?)  =  $  g(s)  ds  = 


1/3 


f  (1  —  e  v  )  +  rjgir}).  Now,  the  stream  function  may  be  written  as 


rv  ,  (-T)2  rY 

Surface  =  /  udj/  =  -  G{rj)  =  /  U  dY’ 

JO  £  lx  do 


21 1/3 


(VI.15) 


The  scalings  are  V  =  —ifix,  U  =  u/P,  and  Y  =  Py. 

Next,  consider  the  wall  velocity  component.  The  stream  function  may  be 
written  as 

i/Vai.  =  -  j*v  dx'  =  -  J*  V  dX'  =  -(40FY)1'4  (VI.16) 


where 


v  = 


T]  = 


r5FlV2  m 


2Y\ 

5  F  11/4 


32  Y3 


x 


61 


The  scalings  are  U  —  ipy,  V  =  v/P,  and  X  =  Px.  Here  F  is  found  from  (V.53)  to 
match  the  surface  flow  to  the  wall  jet. 

For  the  similarity  component  to  the  core  velocity,  ?7sim,  which  is  due  to  the 
wall  jet,  the  stream  function  may  be  written  as 


ipsim  =  ~  [40F]1/4  [n(Y  +  iX)1/A  -(1  +  V2)  S(r  +  iX)1/4] ,  (VI.  17) 


with  the  two  boundary  conditions  that  ip(0,  Y)  =  —  [40F] V1/4  and  ip(X,  0)  =  0. 

For  the  potential  component  to  the  core  velocity  using  the  Green’s  function, 
U green,  the  stream  function  may  be  written  as 


^een  ~Io  2 

where  the  Green’s  function  is 


[(-r)2l 

1/3  dG(X,X o) 

[e2Tx  J 

dy 

y=0 


dx, 


(VI.  18) 


G(X,  Xo)  =  ^~;{ln  [(*  “  Xq)>2  +(y~  y<>)2\  - ln  [(x  “  xo)2  +  (y  +  yo)2] 
-In  [(x  +  x0)2  +  (y  -  yo)2]  +  In  [(x  +  x0)2  +  (y  +  y0)2]  | 


-  so)2  +  (y  -  yo)2]  [(g  +  ^o)2  +  (y  +  y0)2] 
o)2  +  (y  +  yo)2]  [(«  +  ^o)2  +  (y  -  yo)2] 


and 


dGjXjC o), 
dy 


•  _  1 

yo 

yo 

*y=o  7 r 

(x  +  x0)2  +  yl 

(x  —  x0)2  +  yo. 

(VI.  19) 


(VI.20) 


Finally,  summing  all  of  the  components  yields  the  uniform  velocity,  which  can 
be  written  as 


^uniform  (X ,  Y)  =  1psim(X,  Y)  +  V'green(X,  Y)  +  ^surface  (X,  Y/P )  +  ^wall  {X/P,Y) 

'0surface(X,  Oo)  Ipwal l(00>^r)  (VI. 21) 

The  Prandtl  number  P  is  simply  a  parameter  in  the  model  that  represents  how 
thin  the  boundary  layers  are.  The  boundary-layer  model  should  be  independent  of  P. 
Using  the  uniform  stream  function,  a  streamline  field  can  be  plotted.  This  is  shown 
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in  Figure  13.  The  horizontal  velocity  component  increases  as  the  distance  to  the 
wall  decreases,  and  the  vertical  wall  jet  can  clearly  be  seen.  The  flow  into  the  corner 
along  the  surface  turns  down  the  wall,  due  to  the  jet.  The  upper  left  figure  shows  the 
complete  contour  plot  of  the  stream  function  in  the  domain.  The  upper  right  figure 
shows  the  core  stream  function,  Y).  The  middle  set  of  figures  represent  that 

portion  of  the  stream  function  which  is  due  to  the  wall  components,  ^WSLn(X,  Y)  + 
T)— i/^waii (oo,  Y)  on  the  left  and  the  similarity  solution  Y)  on  the  right. 

The  bottom  row  of  figures  shows  those  portions  of  the  uniform  stream  function  which 
are  due  to  the  surface  components  Y)  +  Y/P)  -  i>sw{ace(X,  oo)) 

on  the  left  and  that  due  to  the  Green’s  function  (?/'green(X,  T)  on  the  right.  A  closer 
detail  of  the  corner  stream  function  dynamics  is  shown  in  Figure  14. 

As  can  be  seen,  the  velocity  field  appears  to  be  strong  enough  to  compress 
the  thermal  field  into  the  cold  corner.  The  wall  jet  similarity  solution  dominates 
the  Green’s  function  solution  in  the  core  flow  region  and  overall  in  the  steady-state 
solution.  The  few  streamlines  near  the  corner  in  the  Green’s  function  contour  plot 
suggest  that  h(x)  =  ^  reaches  a  local  maximum  along  the  surface,  not  far 

from  the  corner,  and  the  impact  of  h(x)  is  strong  enough  to  alter  the  velocity  field  in 
this  view.  However,  overall,  h(x)  does  not  appear  to  change  the  expected  behavior  of 
the  velocity  field. 


E.  UNIFORM  TEMPERATURE 

In  a  similar  fashion  to  finding  the  uniform  velocity  approximation,  the  uniform 
temperature  for  the  entire  domain  can  be  written  as 


*  uniform 


(VL22) 


where  Tjnner  is  the  temperature  solution  pertaining  to  the  surface  boundary  layer, 
Touter  is  the  numerically-obtained  temperature  solution  from  the  ADI  method  (solving 
the  unsteady  two-dimensional  heat  equation  with  periodic  velocity  updates),  and 
^match  is  the  matching  condition  which  satisfies  the  flux  condition  at  the  surface. 
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From  Surface 


Steady-Stale  Uniform  Stream  Function 


Steady-State  Uniform  Stream  Function  Due  to  Core  Components 


Figure  13.  Uniform  Stream  Function  and  its  Components 
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Figure  14.  Uniform  Stream  Function  and  its  Components  —  Closer  Detail 
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Recall  that  the  dominant  balance  argument  gives  uTx  —  tyy  as  the  surface  heat 
equation  (see  Chapter  IV). 

Integrating  the  surface  heat  equation  through  the  boundary  layer  gives 

ty(X,y)-ty(X, 0)  =  F uTxdy'  *Tx{X, 0)  f udy' 

J  0  i/0 

=  Tx(X,0)i>(X,y)  =  ty(X,y)  +  TY(X,0)  (VI.23) 

As  y  — »  oc,  ty  —•  0,  and  the  heat  flux  condition  in  the  core  region  becomes 

Ty(X,0)  =  Tx(X,0)4’(X,oo)  (VI.24) 


Solving  for  ty  yields 

roc 

ty(X,y)=Tx(X,0){1>(X,y)-MX,oc)]  =  -Tx(X,0)  udy'  (VI.25) 

Jy 

Integrating  this  expression  from  y  to  oo  gives 

t(X,  oo)  -  t(X,  y)  =  -Tx(X,  0)  J™  (jT  u  dy ")  dy'  (VI.26) 


or 


t(X, y)  =  Tx(X,  0)  r  r  u  dyj'  dy'  (VI.27) 

Jy  Jyf 

Prom  Timman’s  method,  the  horizontal  component  of  velocity  is  u  =  —  j^~rg)rxj  ^ 
where  rj  =  y/6,  6  =  ,  G(y)  =  $  g{£)  d£,  and  g(y)  =  e~s2  ds  =  ^erfcfa). 


Thus, 


ty  —  Tx  V’oo) 

( -T)Txl2/ 3 


VSF 


roc 

/  erfc(r/)  dy' 

Jr) 


tfl'lGW-GH 

(-T)TX 


1  2/3 


ie  v2  -rjg(y) 


(VI. 28) 


66 


and  the  temperature  t  may  be  written  as 


t  = 


Tx  r  (i>  -  M  dy' 

=  -(=f){:§rreicc(ri")dr,"dv'. 

=  (17)  [(1+2^2)s(7?)-??e"V] 


(VI.29) 


Uniform  T*mp*r«tur»  Solution 


As  can  be  seen  in  Figure  15,  the  uniform  temperature  solution  (dotted  contour 
lines)  matches  the  steady-state  solution  (solid  contour  lines)  far  from  the  surface  (as 
77  gets  large),  but  there  is  a  slight  difference  when  77  is  small  (close  to  the  surface), 
due  to  the  G(0)  term.  The  uniform  temperature  prediction  bends  toward  the  surface 
to  intersect  it  perpendicularly,  satisfying  the  Ty  =  0  boundary  condition. 


F.  VORTICITY 

The  vorticity  plays  a  major  role  in  the  cold  corner,  turning  the  flow  from  the 
surface  down  the  wall.  Our  model  has  no  vorticity  in  the  core  flow  region;  however, 
the  numerical  data  shows  some  vorticity  there,  so  the  flow  fields  look  different.  This  is 
due  to  the  hypothesis  that  the  vorticity  that  is  created  far  upstream  from  the  corner 
(along  the  surface)  diffuses  into  the  slow  flow  and  gets  transported  to  the  region  near 
the  cold  corner  and  not  into  the  boundary  layer  regions.  Figure  16  shows  a  contour 
plot  of  the  stream  function  (dotted  contour  lines)  and  vorticity  (solid  contour  lines)  in 
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the  cold  corner,  due  to  Canright’s  numerical  data.  Most  of  the  vorticity  is  confined  to 
thin  regions  along  the  surface  and  wall  (the  boundary  layers),  but  some  vorticity  is  in 
the  core  flow  region  (where  the  streamlines  are  sufficiently  far  apart  to  designate  slow 
flow).  The  main  result  is  that  core  vorticity  is  not  essential  to  the  feedback  process, 
based  on  the  comparison  of  the  thermal  fields  of  the  model  and  numerical  results. 

Vorticity  and  Stream  Function  Contour  Plot 


Figure  16.  Vorticity  and  Stream  Function  Contour  Plot  (due  to  Canright) 

To  quantify  Figure  16,  far  from  the  corner  (x  ~  3),  the  surface  value  of  vorticity 
is  O(103)  smaller  than  than  surface  vorticity  value  at  the  corner  ( x  — >  0).  Recall  that 
a;  =  0atr  =  y  =  0.  Even  at  a  distance  of  <5  from  the  corner  (from  both  the  surface 
and  the  wall),  the  vorticity  value  is  O(102)  smaller  than  the  peak  value  in  the  corner. 
Thus,  the  vorticity  outside  the  boundary  layers  is  non-zero  but  insignificant  compared 
to  the  boundary-layer  values  and  is  not  essential  to  our  model. 
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VII. 


CONCLUSIONS  AND  DISCUSSION 


Complex  Happens.  —  David  Cannght 


The  thermocapillary  stress  condition  at  the  free  surface  gives  a  starting  point 
for  the  scalings  and  dominant  balance.  We  have  obtained  relations  for  the  thermal 
length  scale  l,  the  viscous  thickness  <5,  the  surface  velocity  u,  and  the  core  velocity  U 
scales  in  the  cold  corner,  in  terms  of  the  dimensionless  parameters  M  and  P.  These 
scalings  are 


6  ~  M~lP~\ 
l  ~  M-'P-2, 
u  ~  P, 

U  ~  P2. 


These  scalings  are  consistent  with  the  numerical  results  obtained  by  Canright  in  his 
1994  work  [Ref.  1]. 

The  temperature  gradient,  Tx,  on  the  surface  drives  the  strong  surface  flow, 
similar  to  a  surface  jet.  In  the  viscous  surface  boundary  layer,  a  modification  to 
Timman’s  method  was  employed  which  allowed  the  surface  velocity  to  be  modeled  as 

u  =  e~1/3  (— X1)1/3  (rx)1/3erfc(77),  (VII.  1) 


where  e  is  a  constant  and  77  —  y/8.  The  boundary-layer  thickness  8  in  terms  of  the 
surface  temperature  is 


(VII.2) 


The  surface  flow  then  feeds  into  the  corner  and  into  the  wall  flow,  similar  to  a 
wall  jet.  In  the  viscous  wall  boundary  layer,  a  modification  to  Glauert’s  method  was 
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employed  which  allowed  the  wall  velocity  to  be  modeled  as 

/5  F\1/2 

v={w)  m 


(VII. 3) 


where 


'  5  F  \1/4 
,32  Y3)  '' 


(VII. 4) 


and  F,  the  flux  of  exterior  momentum  flux,  is  found  from  matching  the  surface  velocity 


16y/n  ( Tx 


m  ■ 


(VII.5) 


Using  these  expressions  for  the  velocity  components,  the  stream  functions 
were  calculated.  The  strong  flow  in  the  viscous  boundary  layers  induces  a  relatively 
weak  irrotational  flow  in  the  core  region.  Two  approaches  were  employed  to  find  the 
core  velocities:  a  Laplace’s  equation  solver  using  a  Green’s  function,  and  a  complex- 
variables  solution  modeled  after  a  self-similar  plane  wall  jet.  These  velocities  were 
then  substituted  into  the  unsteady  heat  equation  and  an  alternating-direction  implicit 
(ADI)  numerical  scheme  was  employed  to  solve  for  the  temperature.  The  weak  core 
flow  then  contains  a  thermal  gradient  which  was  used  to  determine  the  Tx  on  the 
surface.  Once  found,  the  surface  temperature  gradient  was  used  to  update  the  surface 
velocity  and  the  process  was  iterated  until  steady  state  was  reached. 

When  the  flow  along  the  wall  is  modeled  as  a  wall  jet,  a  modification  to 
Plotkin’s  second-order  approach  was  used  to  match  our  boundary  conditions  in  the 
quarter-plane.  Thus,  the  potential  solution  becomes 


ip(X,  Y)  =  -4  [lZ(Y  +  iX)1/4  -(l  +  V2)  5(Y  +  iX)1^] .  (VII.6) 

The  velocity  components  may  be  found  by  U  =  and  V  =  —  J^. 

It  is  found  that  the  flow  into  the  cold  comer  is  strong  enough  to  contain 
the  thermal  field.  The  isotherms  are  compressed  along  the  wall  after  steady-state  is 
reached.  The  temperature  profile  obtained  from  the  ADI  method  compares  favorably 
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with  Canright’s  numerical  data  from  1994.  Considering  the  simplicity  of  our  model, 
the  results  are  very  encouraging.  The  thermal  field  results  are  self-consistent. 

The  prediction  model  of  the  uniform  velocity  vector  field  is  encouraging,  in 
that  it  matches  the  inner,  boundary-layer  flow  to  the  outer,  core  flow,  and  subtracts 
the  matching  flow  elements.  The  basic  model  applies  for  the  whole  range  of  Marangoni 
and  Prandtl  numbers  in  the  convective-inertial  regime.  The  uniform  flow  picture  will 
change  with  a  changing  Prandtl  number,  reflecting  the  relative  thickness  of  the  viscous 
boundary  layers,  but  the  underlying  model  does  not  change.  In  addition,  because  this 
model  seems  to  capture  the  cold  corner  feedback  process,  these  results  could  be  useful 
as  a  local  solution  in  a  larger  numerical  problem,  without  the  need  to  numerically 
resolve  the  viscous  boundary  layers. 
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VIII.  FUTURE  WORK 


We  ask,  with  the  extension  of  mathematical  knowledge  will  it  not  finally 

BECOME  IMPOSSIBLE  FOR  THE  SINGLE  INVESTIGATOR  TO  EMBRACE  ALL  THE  DEPARTMENTS 
OF  THIS  KNOWLEDGE?  —  David  Hilbert 


There  are  a  few  limitations  to  this  work.  The  grid  used  in  the  numerical 
approximations  is  uniform.  To  get  a  more  accurate  description  of  the  velocity  and 
temperature  fields,  a  non-uniform  grid  should  be  employed,  which  is  finer  near  the 
corner  and  more  coarse  away  from  the  corner,  as  in  Canright’s  1994  work.  Even  so, 
the  current  data  shows  good  correlation  with  the  previous  work.  A  nonuniform  grid 
would  require  a  more  complicated  algorithm  but  may  run  faster  for  a  given  resolution 
in  the  corner,  especially  since  running  of  the  current  ADI  codes  with  updates  required 
hours  for  several  thousand  iterations. 

More  effort  should  be  placed  on  the  effects  of  vorticity  in  the  corner.  In  the 
numerical  data,  the  vorticity  does  not  appear  to  be  entirely  confined  to  the  boundary 
layers,  yet  the  vorticity  in  the  core  flow  region  seems  to  have  no  significant  impact 
on  the  temperature  field  solution.  It  appears  that  the  core  flow  does  not  turn  down 
the  wall  if  vorticity  is  not  present.  The  vorticity  in  the  surface  and  wall  boundary 
layers  must  be  strong  enough  to  force  the  nonlinear  flow  at  the  free  surface  down  the 
rigid  wall.  At  this  point,  we  can  estimate  the  magnitude  of  the  vorticity,  but  more 
research  is  needed  to  determine  how  to  combine  the  vorticity  into  the  potential  flow.  A 
meaningful  approach  to  modeling  the  vorticity  throughout  the  domain  would  provide 
a  good  understanding  of  the  rotation  mechanism  in  the  flow.  This  was  addressed  in 
Chapter  VI. 

This  research  simplified  the  physical  problem.  The  geometry  was  revised,  some 
material  properties  were  ignored,  and  the  free  surface  was  assumed  flat  due  to  the 
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amount  of  surface  tension.  If  a  geometry  change  occurs,  that  should  be  investigated 
in  future  research. 

Also,  higher-order  modifications  to  Timman’s  method  should  be  attempted. 
This  would  provide  more  accuracy  near  the  free  surface.  The  current  model  predicts 
the  decay  of  the  velocity  through  the  boundary  layer  fairly  well,  until  the  interface 
with  the  core  region  is  approached,  and  then  the  Timman  modification  decays  faster. 
If  a  four-parameter  model  could  be  developed,  it  might  better  depict  the  velocity 
profile  using  the  thermocapillary  condition  at  the  surface. 

As  more  experiments  are  conducted  in  the  low-gravity  environment,  such  as 
those  using  the  NASA  space  shuttles  by  Kamotani  and  Ostrach  and  others,  a  valida¬ 
tion  of  the  current  work  using  experimental  techniques  should  be  attempted. 


74 


APPENDIX  A.  TIMMAN’S  METHOD  ALONG 

THE  SURFACE 

Even  though  the  one-parameter  model  is  very  encouraging,  the  goal  is  to  try 
to  derive  a  more  accurate  velocity  profile  than  the  simplest  approach  in  Chapter  V, 
since  the  surface  tension  and  thermocapillary  condition  at  the  surface  is  driving  the 
flow.  Now,  suppose  the  velocity  profile  for  u(y)  (where  the  A  dependence  is  implicit) 
is 

u(y)  =  f(tj)  =  e"*2  (a{X)  +  c{X)t2)  dt  -  b{X)e~Tt\  (A.l) 

This  introduces  a  damping  term  into  the  velocity  profile.  Is  there  a  more  accurate 
expression  for  the  boundary-layer  thickness?  At  the  surface,  y  =  r)  =  0,  which  implies 
that 

ti(0)  =  /(0)  =  e~t2  (a(X)  +  c(X)t2)  dt  -  b(X) 

-  ^«P0  +  ^o(X)  -  b(X)  (A.2) 

In  the  core  flow  region,  as  y  —*■  oo,  77  — >  00,  and  /( 00)  =  0. 

Differentiating  f'(rj )  =  e^2  [-a(X)  -  c(X)r)2  +  2b(X)r]}.  Evaluating 

this  at  the  surface  gives  /'( 0)  =  -a,  which  yields  a(X)  =  -8(X)Tx,  just  as  in 
Chapter  V  (when  6(A)  =  c(X)  =  0). 

Differentiating  further  yields  f"(r)),  which  when  evaluated  at  the  surface  gives 
/"( 0)  =  26(A).  Since  f(rj)  =  u,  f'(rj)  =  8(X)uy,  and  f"{r})  =  82(X)uyy.  Thus, 
6(A)  =  1 62(A)  u( 0)  ttx(0)  at  the  surface. 

Taking  the  third  derivative  of  /(??),  /"'(r?)  is  a  very  lengthy  expression.  uyyy 
can  be  rewritten  as  ~  (uy y),  which  after  substitution  into  the  momentum  equation 
becomes  |(uu,  4-  V  uy).  This  can  be  simplified  using  the  continuity  equation  and 
the  fact  that  r(0)  =  0.  Evaluation  at  the  surface  yields  /"'( 0)  =  <53  u(0)  TXx-  Since 
/"'( 0)  =  2a(X)  —  2c(A),  we  can  solve  for  c(A). 


75 


The  proposed  velocity  profile  at  the  surface  is  therefore 

«(»)  =  /(!?)  =  -  r  [t(X)Tx  +  { STx  +  i«3u(0)Txx}  f2]  e-‘2dt-^S1(X)u(0)ux(0)e-’\ 

(A-3) 

Hence,  a  first-order,  nonlinear,  nonhomogeneous  ordinary  differential  equation  for  the 
velocity  can  be  developed  in  terms  of  the  boundary-layer  thickness  and  the  tempera¬ 
ture  gradient  as  (let  u(0)  =  uq): 

«o«o  +  “0  (I  +  X57ix)  =  "■ 2TTx  (A'4) 

where  u'0  =  ^u(O).  This  gives  one  equation  with  two  unknowns,  uq(X)  and  S(X). 

We  also  have  <52  =  —T  =  6  /0°°  u2d  rj.  This  type  of  ordinary  differential  equation  (A.4) 
can  be  known  as  an  Abel  Equation  of  the  Second  Kind  [Ref.  33].  Finding  a  solution 
to  this  type  of  equation  is  available  only  for  certain  cases,  such  as  if  the  coefficients  of 
uqu'q  and  uo  fall  into  special  categories  outlined  in  [Ref.  34,  35].  Unfortunately,  this 
particular  ODE  does  not  fall  into  one  of  those  special  cases. 

A  second  approach  to  modeling  the  velocity  with  a  Timman-type  method  of 
two  parameters  is  algebraic.  Write  the  velocity  profile  as 

u  =  f(rj)  =  f  ae~t2dt  —  be*7? . 

JT) 

(Remember  that  the  variables  a  and  b  are  functions  of  X;  however,  for  simplicity, 
the  a(X)  and  b(X)  will  be  replaced  by  a  and  b ).  Using  three  derivatives  evaluated 
at  y  =  0  will  enable  us  to  obtain  an  expression  for  u  in  terms  of  the  boundary-layer 
thickness  and  then  in  only  temperature  terms.  Recall  that 

m  =  -  b 

Taking  derivatives  of  /  with  respect  to  77  and  using  the  chain  rule  since  y  —  6tj,  we 
obtain 

f(y)  =  e~v2(-a  +  2bq)  =  8uy 

f(0)  =  —a  =  6  (uy)y=0  =  STx  (A.5) 
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Thus,  a  =  — 6Tx •  Continuing, 


f"(v)  =  e~rt\2b  +  2ar]  -  Abrf)  =  82uyy 
/"( 0)  =  2b 

f'(v)  =  e-T>2[2a-8brj-2r](2b  +  2aV-4bn2)}=63uyyy 
0)  =  2a  —  —26  Tx 

Recall  that  uyyy  =  ^  (uyy)  =  (uux  +  Vuy )  =  u  Txx ■  Thus 

/'"( 0)  =  63u(0)Txx  =  63  (£.  -  =  «3  (-^Tx  -  b\  Txx 

which,  when  equating  the  /'"( 0)  terms,  finally  yields 

=  +  (A.8) 


(A.6) 

(A.7) 


so  that 


6  = 


2Tx  _  V^Lcrp 
62TXX  2 


(A.9) 


and  the  velocity  profile  can  now  be  modeled  as 

2  Tx 


u 


=  -  8  Txe~t2dt- 

J  n 


V* cT 

6*TXX~  TSTx 


o-*r 


8  Tx  erfc(7/) 


2  Tx  v^F 


82T; 


8T- 


X 


XX 


o-r 


Recalhng  (V.9),  the  boundary-layer  thickness  8  can  be  solved  as 

«(*)  =  t 

Jo  u  dr! 

This  yields  a  quadratic  equation  in  <53,  which  is: 


86  \T\ 


+  s3 


(A.  10) 


(A.11) 


r-S(vM+3l?-0  <A-12> 

At  the  corner,  with  T,  Tx,  and  Txx  values  taken  from  the  numerical  data 
at  x  =  0,  the  boundary-layer  thickness  8(X)  is  0(M~1P~ 1),  as  expected.  As  X 
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increases  from  the  wall,  6(X)  decreases  in  value,  also  as  expected.  However,  after 
a  certain  value  of  X,  still  within  the  comer  region,  the  solution  to  (A.12)  results  in 
complex  roots.  Unfortunately,  this  means  that  the  algebraic  two-parameter  approach 
had  to  be  abandoned. 
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APPENDIX  B.  CALCULATION  OF  e 


In  the  surface  boundary-layer  thickness  equation,  (V.22),  the  constant 

7T 

£=4Jo  terfc^)]  ^ 

is  found  to  be 

•-t('-t) 

Here  is  the  work.  Integration  by  parts  is  necessary.  Define  I  =  /0°°  [erfc(r/)]2  dr],  such 
that  e  =  jl.  Let  u  =  [erfc^)]2  and  dv  =  dr].  Then  du  —  erfc(r?)  e-7?2)  dr]  and 


v  —  7],  Thus 


Since  erfc(oo)  =  0, 


I  —  rj  [erfc(77)]' 


00  4  f°°  __  2 

\ — 7=  /  Ve  77  &k(r))dr) 

7T  JO 


0 


4  roo  2 

;=  /  rje  77  ei{c(r))dr) 
7T  Jo 


2  r°°  2 

=  —7=  /  2rje  71  eiicMdr) 

y/TT  JO 

Integrate  by  parts  again.  Let  u  =  erfc(77)  and  dv  =  2r]  e~^  dr\.  Then  du  =  — e~n2  dr] 


and  v  =  —e  n  .  Hence, 


1  = 


2 


-2-Te^dr] 

y/l r  Jo 

Use  the  transformation  tu  =  \/2r/.  Then  dw  —  \/2d?].  Substituting, 

2  1  r°°  2 

e~w  dw 

\TK  y/2  JO 


I  = 


V* 

2 

V5F 


'  V5  V*]  2  /  _  V2\ 

2  J  0F  l  2  ) 
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Therefore, 


f(-f) 
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APPENDIX  C.  VORTICITY  FLUX 


The  vorticity  equation  in  the  wall  boundary  layer  is 


U  Ux  +  V  OJy  =  Uxx 


This  can  be  rewritten  as  (U  ui)x  +  (vu)y  =  wxx  using  continuity.  Now,  integrate  this 
with  respect  to  x,  obtaining 

d  f°° 


o°  d  r°° 

U  u>  +  /  vujdx  =  ua 

x=o  dY  Jo 


oo 

x=0 


(C.l) 


As  x  — ►  oo,  u  —*■  0  (and  ux  —*■  0),  and  when  x  =  0,  U  =  0,  so  that 

d  r°° 


d  r00 

——  /  vudx  =  —ojx 

dY  Jo  x=° 


(C.2) 


The  vorticity  can  also  be  written  as  u  =  vx—Uy -  Prom  the  continuity  equation, 


U  =  - 


Vy  dx 


which  yields 

—  Vx  +  —  Jq  Vy  dx  =  Vx  +  Jq  Vyy  dx.  (C.3) 

A  scahng  analysis  is  performed  on  the  vorticity  terms.  In  the  wall  boundary 
layer,  x  ~  A  ~  x,  y  ~  L  ~  M~lP~2Y ,  u  ~  P2  U,  v  ~  Pv,  and  u  ~  MP2 u. 

Substituting  into  u  —  vx  —  Uy  yields  MP2u>  =  MP2  vx  +  MP i  Uy.  The  last  term 
on  the  right  is  much  smaller  than  the  first  term  on  the  right  (since  Pel),  giving 
u  =  vx.  Thus,  the  vorticity  u  can  be  characterized  by  vx,  ignoring  the  /q  vyy  dx  term 
as  small.  This  result  is  consistent  with  the  numerical  data  of  Canright  [Ref.  1],  as 
shown  in  Figure  17  (the  vorticity  plot  is  the  solid  line,  the  velocity  derivative  plot 
is  the  dotted  line).  As  can  be  seen,  the  two  plots  are  virtually  identical  for  given 
horizontal  slices  of  the  vorticity  and  vertical  velocity  flux  (in  this  case,  at  a  depth  of 
U  =  0.0584). 
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Vorticity  vs  Velocity  Flux 


Figure  17.  Typical  Vorticity  and  Vertical  Velocity  Flux  Profiles  at  the  same  Depth 
from  the  Free  Surface 

Back  to  the  equation  for  the  vorticity  flux, 


=  0  (C.4) 

This  shows  that  the  vorticity  flux  at  the  wall  is  zero!  Very  unexpected.  Sub¬ 
stituting  in  the  expression  for  v  and  vx  and  then  evaluating  the  integral  does  in¬ 
deed  yield  0.  Perhaps  a  higher  order  expansion  for  the  velocity  is  needed,  such  as 
v  «  Pvi  +  P2v2  - ,  to  resolve  the  diffusion  of  vorticity  from  the  wall. 
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APPENDIX  D.  TIMMAN’S  METHOD  ALONG 

THE  WALL 


The  velocity  profile  approximation  for  the  wall  boundary  layer  is  a  little  more 
complicated  than  the  surface  profile  approximation.  The  no-slip  condition  at  the  wall 
requires  that  the  velocity  is  zero  when  x  =  0.  Also,  the  core  velocity  goes  to  0,  so 
the  boundary  condition  as  x  — *  oo  requires  a  zero  velocity.  Between  these  two  values 
of  x,  the  velocity  increases  to  the  point  in  the  boundary  layer  where  the  positive 
vorticity  goes  to  zero,  then  the  velocity  decays  through  the  boundary  layer  where  the 
vorticity  is  negative  (that  vorticity  which  is  convected  down  from  the  surface)  until 
the  velocity  decays  to  zero. 

The  following  approximate  velocity  profile  is  assumed: 

v(x)  =  f(r})  =  -  f  e~t2(a  4-  ct 2)  dt  +  e-7?2(6  +  drj2),  (D.l) 

Jr) 

where  f(rj)  =  v,  the  vertical  (dominant)  component  of  the  wall  velocity,  and  77  = 
where  A  is  the  thickness  of  the  wall  boundary  layer.  The  boundary  conditions  are: 
at  x  =  0, 77  =  0,  and  f(r])  =  0,  due  to  the  no-slip  condition.  As  x  — >  00,  77  — >  00,  and 
/(» l)  0- 

The  determination  of  the  coefficients  a,  b,  c,  and  d  must  be  made.  All  four 
coefficients  are  functions  of  Y,  the  vertical  direction.  Start  with  the  momentum 
direction  in  the  Y-direction, 

Uvx  +  VVy  =  vxx  (D.2) 

Integrate  this  equation  with  respect  to  x,  from  x  =  0  at  the  wall  to  x  =  00,  outside 
the  vertical  boundary  layer.  This  yields 


r°o  I 00 

/  [Uvx  +  vvY)dx  =  vx\ 

Jo  10 

Substitute  the  shear  stress  at  the  wall,  r0,  as  r0  =  =  vx.  Thus, 

roc 

—  To  =  /  ( Uvx  +  Wy)dx 

JO 


(D.3) 


(D.4) 
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(Uv) 


+ 


dx 


d  rcc 

dY 


r  oo 

/  »= 

Jo 


dx 


Define  the  momentum  thickness  to  be  A2,  where 

fCC 

A2  =  /  v 2  dx. 

J  x=0 


Then 


roc  roo 

A2  =  /  v2dx  =  f2  dx. 

Jx— 0  Jx=0 


Since  77  =  4,  dx  =  A  dr/,  and 


/•oo 

A2  =  A  /  v2  dr], 

Jx—O 


with  A  being  the  boundary-layer  thickness.  Therefore, 

dA2 


-t0  = 


(D.5) 


(D.6) 


(D.7) 


(D.8) 


The  derivative  conditions  on  the  function  f(rj)  evaluated  at  the  wall  are  as 
follows: 

'<">  -  s-x(“+ i) 

/'(O)  =  a 

/"( 0)  =  2(d  —  6) 

.HO)  =  2(c  —  a) 

/iv(0)  =  12(6 -2d) 


Since  /  =  v,  f(0)  =  0,  as  the  velocity  at  the  wall  is  zero.  Since  there  is  no  pressure 
gradient  term  in  the  momentum  equation,  the  pressure  is  constant,  and  /"( 0)  =  0. 
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In  addition,  /"'( 0)  =  0,  from  the  given  boundary  conditions  [Ref.  31].  These  last 
two  conditions  yield  a  =  c  and  b  =  d.  The  first  condition  gives  a  =  3^b-  0ne 
more  condition  is  needed  to  solve  for  the  four  coefficients.  The  value  of  the  shear 
stress  at  the  wall  is  unknown,  so  /'( 0)  is  not  specified,  and  the  expression  for  /iv( 0)  is 
needed.  Taking  the  fourth  derivative  of  /  with  respect  to  rj  yields  f  iv(ri )  =  A4vxvxy. 
Since  f'(rj)  =  Avx,  then  a  =  Avx,  evaluated  at  x  =  0,  and  the  last  condition  yields 
/*<»)-  A3o^(f). 

Equating  the  various  derivative  conditions  yields  two  nonlinear  equations  in¬ 
volving  the  boundary-layer  thickness  A  and  the  coefficient  a.  These  are 


-m=ASaw(i) 

Evaluating  the  integral  associated  with  A2  numerically,  the  first  equation  can  be 
worked  to  the  form 


-To  = 


-Vx 


a 

A 


d 

dY 


(r2A) 


—  Ay  +  2aAay'j 


£ 

A 


—  aAy  -f-  2Afl y 


while  the  second  equation  can  be  worked  to  a  similar  form 


-126  =  — = 


(D.9) 
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9  Vtt 

"A" 


aAy  —  A<2y 


(D.10) 


These  equations  may  be  further  reduced  to 

ay  =  “  3v^) 

(D-ll) 

and 

Ay  =  (1  +  6  v/tt) 

Z\a 

(D.12) 

Initial  conditions  need  to  be  specified  on  the  displacement  thickness  (call  it 
the  mass  flux)  and  on  the  flux  of  vorticity.  The  mass  flux  and  the  vorticity  flux  must 
be  matched  to  the  free  surface  boundary  layer,  so  more  than  one  free  parameter  may 
be  needed. 

If  d  =  0,  then  from  the  initial  conditions,  6  =  0,  which  further  means  that 
a  =  —  f .  The  velocity  profile  is  then 

ro o  0 

/(? 7)  =  v  =  a  /  e_t  (2 12  —  1)  dtt.  (D.13) 

Jr\ 

This  approach  proved  to  be  more  difficult  than  expected  and  was  abandoned 
in  favor  of  the  modification  to  Glauert’s  method. 
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APPENDIX  E.  NUMERICAL  CODES 


.  May  YOU  CONVERGE  WITHOUT  DELAY.  -  Dan  Zinder 


The  computer  programs  listed  here  axe  supplied  on  an  “as  is”  basis,  with  no 
warrantees  of  any  kind.  The  author  bears  no  responsibility  for  any  consequences  of 
using  this  program. 
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/*  this  file:  fluid. h  -  header  definitions  for  cold  corner  */ 

#include  <stdio.h> 

#include  <math.h> 

#define  maxx  3001 
#define  maxxvec  (maxx  -  2) 

#define  Ma  10000 
#define  Pr  0.01 

#define  eps  (sqrt (3. 14159265) *(1.0-sqrt (2. 0)/2.0)/2.0) 

#ifdef  MAIN 

int  i,  j,  k,  nx,  ny; 
int  size; 

double  hx,  hy,  dt; 

double  x [maxx] ,  y [maxx] ,  b [maxx] ,  xx [maxx] ,  f  [maxx]  ; 
double  xconl,  xcon2,  yconl,  ycon2; 

double  cl  [maxx]  [maxx] ,  c2,  c3[maxx]  [maxx]  ,  c4[maxx]  [maxx]  ,  c5, 
c6[maxx]  [maxx]  ; 

double  dl  [maxx]  [maxx]  ,  d2 ,  d3  [maxx]  [maxx]  ,  d4  [maxx]  [maxx]  ,  d5 , 
d6 [maxx] [maxx] ; 

double  t  [maxx]  [maxx]  ,  temp  [maxx]  [maxx]  ,  tstar  [maxx]  [maxx]  ; 

double  hofx[maxx],  integrand  [maxx]  ,  gmfun[maxx]; 

double  gpsi[maxx]  [maxx]  ,  jpsi[maxx]  [maxx]  ,  psi[maxx]  [maxx]  ; 

double  wpsi  [maxx]  [maxx]  ,  spsi  [maxx]  [maxx]  ,  unifpsi  [maxx]  [maxx]  ; 

double  u [maxx] [maxx] ,  v [maxx] [maxx] ; 

double  surfu,  massflux,  F,  eta,  g; 

double  delta [maxx] ,  seta [maxx] ; 

double  tempa[maxxvec]  ,  tempb  [maxxvec]  ,  tempo  [maxxvec]  ; 
double  tstarn[maxx]  ,  Tx[maxx],  Txx[maxx]  ; 
double  A [3] [maxx] ,  C  [3] [maxx] ; 

#undef  MAIN 
#else 

extern  int  i,  j,  k,  nx,  ny; 

extern  int  size; 

extern  double  hx,  hy,  dt; 

extern  double  x  []  ,  y  []  ,  b  []  ,  xx  []  ,  f  []  ; 

extern  double  xconl,  xcon2,  yconl,  ycon2; 

extern  double  cl  []  [maxx],  c2,  c3  []  [maxx]  ,  c4[][maxx],  c5,  c6[]  [maixx]; 
extern  double  dl  []  [maxx]  ,  d2 ,  d3  []  [maxx]  ,  d4  []  [maxx]  ,  d5 ,  d6  []  [maxx]  ; 
extern  double  t[][maxx],  temp  []  [maxx]  ,  tstar  []  [maxx]  ; 
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extern  double  hofx[]  ,  integrand  []  ,  grnfun[]  ; 

extern  double  gpsi  []  [maxx]  ,  jpsi  []  [maxx]  ,  psi  []  [maxx]  ; 

extern  double  wpsi  []  [maxx]  ,  spsi  []  [maxx] ,  unifpsi  []  [maxx] ; 

extern  double  u  []  [maxx]  ,  v  []  [maxx]  ; 

extern  double  surfu,  massflux,  F,  eta,  g; 

extern  double  delta  []  ,  seta[]  ; 

extern  double  tempa[]  ,  tempb[],  tempo  []; 

extern  double  tstam[]  ,  Tx[]  ,  Txx[]; 

extern  double  A  []  [maxx]  ,  C  []  [maxx]  ; 

#endif 
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/*  this  file:  main.c  -  main  program  for  cold  corner  */ 

#define  MAIN 
#include  "fluid. h" 

#define  LINESIZE  128 

#define  eps  (sqrt (3. 14159265) *(1 .0-sqrt (2.0) /2.0)/2 .0) 

main (argc , argv) 

int  argc ;  char  *argv [] ; 

•C 

double  adi(),  tol=0.001,  norm; 
int  n,  m,  ninner=6,  nouter=6; 

if  (  argc  >  1  )  sscanf(  argv[l]  ,  "‘/.d",  feninner  ); 

if  (  argc  >  2  )  sscanf  (  argv  [2]  ,  "’/.d" ,  ftnouter  ) ; 

if  (  argc  >  3  )  sscanf  (  argv [3] ,  "yilf",  &tol  ); 

n  =  10; 
adiinit () ; 

for  (m  =  1;  n  >  2  &&  m  <=  nouter;  m++){ 
tx() ; 

strength  () ; 
jet() ; 
green () ; 

for  (i  =  0;  i  <=  nx;  i++){ 
for  (j  =  0;  j  <=  ny;  j++){ 

psi  [i]  [j]  =  jpsi  [i]  [j]  +  gpsi  [i]  [j]  ; 

> 

> 

for  (i  =  0;  i  <=  nx;  i++){ 
for  (j  =  0;  j  <=  ny;  j++M 
if  (i  ==  0  &  j  !=  0){ 

u[i]  [j]  =  Cpsi [i]  [j+1]  -  psi [i]  [j-1])  /  C(2.0)*(y[j+1]  -y[j])); 

v[i]  [j]  =  -  (psi[i+l]  [j]  -  psi [i]  [j] )  /  (x[i+l]  -  x[i]);> 

else  if  (j  ==  0  &  i  !=  0){ 

u[i]  [j]  =  (psi  [i]  [j+1]  -  psi  [i]  [j] )  /  (y[j+l]  -  y[j]); 
v[i][j]  =  -  (psi[i+l][j]  -  psi[i-l][j])  /  C(2.0)*(x[i+1]  -  x[i]));> 
else  if  (i  ==  nx){ 

u[i]  [j]  =  (psi [i]  [j+1]  -  psi [i]  [j-1] )  /  (2.0*(y[j+l]  -  y[j])); 

v[i]  [j]  =  -  (psi  [i]  [j]  -  psi  [i-1]  [j] )  /  (x[i]  -  x[i-l]);> 
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else  if  (j  ==  ny){ 

u[i][j]  =  (psi  [i]  [j]  -  psi  [i]  [j-1] )  /  (y[j]  -  y[j-l]); 

v[i]  [j]  =  -  (psi  [i+1]  [j]  -  psi  [i-1]  [j] )  /  (2.0*(x[i+l]  -  x[i])); 

u  [nx]  [ny]  =  v  [nx]  [ny]  =  0 . 0 ;  > 

else{ 

u[i][j]  =  (psi [i]  [j+1]  -  psi [i]  [j-1])  /  ((2.0)*(y[j+l]  -  y[j])); 
v[i]  [j]  =  -  (psi  [i+1]  [j]  -  psi  [i-1]  [j])  /  ((2.0)*(x[i+l]  -  x[i]));> 
/*  printf  (  "7.1. 2f  7.1. 2f  %1.5f  7.1.5f\n",  x[i],  y[j],  u[i][j],  v[i][j]>;  */ 

> 

> 

norm  =  1.0; 

for  (n  =  1;  norm  >  tol  &&  n  <  ninner;  n++){ 
norm  =  adi() ; 

/*  printf  ("Iteration  number  7.d,  norm:  7.f  \n"  ,n,norm)  ;  */ 

> 

f printf  (stderr,  "Iteration  =  7.d\n",  m)  ; 

> 

for  (i  =  0;  i  <=  nx;  i++){ 
for  (j  =  0;  j  <=  ny;  j++){ 

printf  (" 74. 3f  %1.3f  7.1.6f\n",  x[i]  ,  y[j],  temp[i][j]); 

> 

> 
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/*  this  file:  adiinit.c  -  initialization  for  cold  corner  */ 
#include  "fluid. h" 
adiinitO  { 


nx  =  ny  =  maxx  -  1; 
hx  =  hy  =  3.0  /  nx; 
dt  =  0.4  *  hx; 

xconl  =  dt  /  (hx  *  hx) ; 
xcon2  =  dt  /  (4 . 0  *  hx) ; 
yconl  =  dt  /  (hy  *  hy) ; 
ycon2  =  dt  /  (4.0*  hy) ; 

/*  printf  ("%1.4f  71. 4f  7,1.4  %1.4f\n",  xconl,  xcon2,  yconl,  ycon2) ;  */ 

c2  =  1.0  +  xconl; 

c5  =  1.0  -  yconl ; 

d2  =  1.0  +  yconl ; 

d5  =  1.0  -  xconl ; 

for  (i  =  0;  i  <=  nx  ;  i++){ 
x[i]  =  y[i]  =  i  *  hx; 
f[i]  =  0.0; 

> 

/*  printf ("  initial\n"); 

printf  ("x[i]  y[j]  temp[i]  [j]\n")  ;  */ 
for  (i  =  0;  i  <=  nx  ;  i++){ 
for  (j  =  0;  j  <=  ny  ;  j++){ 
if  (i  ==  0) 

temp[i]  [j]  =  tstar[i]  [j]  =  -1.0; 
else  if  (i  ==  nx) 

temp[i]  [j]  =  tstar[i][j]  =  0.0; 
else  if  (j  ==  0) 

temp[i]  [j]  =  -exp(-10.0  *  x[i]); 
else 

temp[i][j]  =  tstar[i][j]  =  0.0; 

/*  printf  ("7,1. 3f  7,1. 3f  7,1.4f\n",  x[i],  y[j]  ,  temp[i][j]);  */ 

> 

> 

} 
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/*  this  file:adiinit2.c  -  initialization  for  adi  solver  at  fine  grid.  */ 


#include  "fluid. h" 
#def ine  LINESIZE  128 

adiinitO  { 


nx  =  ny  =  maxx  -  1; 
hx  =  hy  =  3.0  /  nx; 
dt  =  0.4  *  hx; 

xconl  =  dt  /  (hx  *  hx) ; 
xcon2  =  dt  /  (4.0*  hx) ; 
yconl  =  dt  /  (hy  *  hy) ; 
ycon2  =  dt  /  (4.0  *  hy) ; 

/*  printf  ("‘/,1.4f  7,1. 4f  7.1.4  7.1.4f\n",  xconl,  xcon2,  yconl,  ycon2) ;  */ 

c2  =  1.0  +  xconl; 
c5  =  1.0  -  yconl ; 
d2  =  1.0  +  yconl ; 
d5  =  1.0  -  xconl; 

for  (i  =  0;  i  <=  nx  ;  i++){ 
x[i]  =  y[i]  =  i  *  hx; 
f [i]  =  0.0; 

} 

for  (i  =  0;  i  <=  nx;  i++){ 
for (j  =  0;  j  <=  ny;  j++){ 

scanf  ("%*lf 7.*lf7.1f "  ,  &(temp[i]  [j]))  ;  /*  Read  in  converged  data  */ 

tstar[i][j]  =temp[i][j]; 

> 

> 

for  (j  =  0;  j  <=  ny  ;  j++){ 

temp[0]  [j]  =  tstar[0]  [j]  =  -1.0; 
temp[nx]  [j]  =  tstar[nx]  [j]  =  0.0; 

> 

> 
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/*  this  file:  tx.c  This  file  calculates  the  thermal  gradients, 

Tx  and  Txx,  and  the  heat  flux.  */ 

#include  "fluid. h" 

tx() 

for  (i  =  0;  i  <=  nx;  i++M 
if  (i  ==  0) 

Tx[i]  =  (temp  [i+1]  [0] -temp  [i]  [0])  /  (x[i+l]-x[i]) ; 
else  if  (i  ==  nx) 

Tx[i]  =  (temp  [i]  [0] -temp  [i-1]  [0])  /  (x[i]-x[i-l])  ; 
else 

Tx[i]  =  (temp [i+1]  [0] -temp [i-1]  [0]  )  /  (2.0*  (x[i+l]-x[i]>)  ; 
/*  printf  ("7,1.2f  %1.4f  7.1.4f\n",  x[i]  ,  temp[i][0],  Tx[i]);  */ 

> 

for  (i  =  0;  i  <=  nx;  i++){ 
if  (i  ==  0) 

Txx  [i]  =  (Tx[i+1]  -  Tx[i]>  /  (x  [i+1] -x  [i]  ) ; 
else  if  (i  ==  nx) 

Txx[i]  =0; 
else 

Txx[i]  =  (temp  [i+1]  [0] -2*temp[i]  [0]+temp  [i-1]  [0]  )  / 
pow((x[i+l]-x[i]),2.0); 

/*  printf  ("%1.2f  %1.4f  /Cl.4f\n",  x[i]  ,  temp[i][0],  Txx[i]);  */ 

> 

for  (i  =  0;  i  <=  nx;  i++){ 

f[i]  =  (-1. 0/2.0)*  pow((fabs (-temp [i]  [0] *Tx[i] )/eps) , (2. 0/3.0) ) ; 

> 
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/*  this  file:  strength2.c  This  is  a  file  to  determine  the 
characteristic  velocity  from  the  surface  boundary  layer  into  the  corner.  */ 

#include  "fluid. h" 

strength () 

{ 

F  =  f abs ( ( 16 . 0/5 . 0) * (sqrt (3 . 14159265) ) * (pow(Tx [0] / (eps*2 . 0)  ,  1 . 0/3 . 0) ) ) ; 

/*  printf  ("F  =  7.1.6f\n",  F) ;  */ 

> 
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/*  this  file:  jet.c  This  file  determines  the  velocity  from  the 
potential  flow  due  to  a  wall  jet  in  the  quarter-plane .  */ 

#include  "fluid. h" 

jetO 

double  theta,  con; 
con  =  (5. 0/2.0)  *  F; 

/*  printf("x[i]  yCj]  jpsi  [i]  [j]\n")  ;  */ 
for  (i  =  0;  i  <=  nx;  i++){ 
for(j  =  0;  j  <=  ny;  j++){ 
theta  =  atan2(x[i]  ,y  [j]  )  ; 
if  (j  ==  0) 

jpsi[i]  [j]  =0.0; 
else  if  (i  ==  0) 

jpsi [i] [j]  =  -2.0  *  pow(con,  0.25)  *  pow(y [j] ,0.25) ; 
else 

jpsi [i]  [j]  =  -2.0  *  pow(con,  0.25)  *  pow((x[i]*x[i]+y[j]*y[j]) , 
0.125)  *  (cos (thet a/4.0)  -  (1 . 0+sqrt (2. 0))*s in (theta/4.0)) ; 

/*  printf  ('7,1. 2f  7.1. 2f  7,1.4f\n",  x[i],  y[j]  ,  jpsi [i]  [j] )  ;  */ 

> 

> 

} 
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/*  this  file:  green. c  This  file  solves  Laplace’s  equation  in  core.  */ 


#include  "fluid. h" 

greenO 

{ 

double  piece; 

for  (i  =  0;  i  <=  nx;  i++){ 
if  (i  <  nx) 

hof x  [i]  =  (  pow(  ( (-temp [i]  [0] ) * (-temp [i] [0]  ) / (eps*eps*Tx [i] ) ) , 

1. 0/3.0)  )  /  (-2.0); 
else 

hofx[i]  =0.0; 

/*  printf  ("7,1. 2f  */,1.12f\n",  x[i],  hofx[i]);  */ 

} 

for  (i  =  1;  i  <=  nx;  i++){ 
for(j  =  1;  j  <=  ny;  j++){ 
piece  =  0.0; 
integrand [0]  =  0.0; 
for  (k  =  1;  k  <=  nx;  k++){ 

gmfun[k]  =  (  y[j]/((x[k]+x[i])*(x[k]+x[i])  +  y[j]*y[j])  - 
y C j ] / ( (x [k] -x [i] ) * (x [k] -x [i] )  +  y[j]*y[j])  )/(3. 14159265)  ; 
integrand  [k]  =  (-  grnfunCk]  *  hofx[k]); 

piece  +=  C(x[k]-x[k-l]  )  /  2.0)  *  (integrand [k]  +integrand [k— 1]  )  ; 

/*  printf  ("7,1. 4f\n",  integrand  [k]  )  ;  */ 

> 

gpsi  [i] tj]  =  piece  +  (hofx[0]  *  (2.0/3.14159265)  *  atan2(y [j] ,x[i])) ; 

> 

> 

for  (i  =  0;  i  <=  nx;  i++){ 
gpsi[0]  [i]  =  hofx[0]  ; 
gpsi[i]  [0]  =  hofx[i]  ; 

} 

> 
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/*  this  file:  adi.c  -  solve  the  heat  equation  for  the  potential  flow 
due  to  a  wall  jet  in  the  cold  corner  (quarter  plane). 

2  2 


d  T 

d  T 

d  T 

d  T 

d  T 

M( -  + 

u -  + 

v - ) 

= -  + 

— 

d  t 

d  x 

d  y 

2 

2 

d  x  d  y 

with  BC: 

at  y  =  0:  Ty  =  0,  uy  =  Tx,  v  =  0 

at  x  =  0:  T  =  -1,  u  =  v  =  0 

using  the  symmetric  ADI  method,  which  uses  two  subprograms  for  computing 
the  L  U  decomposition  of  a  tridiagonal  matrix  and  then  for  solving 
L  U  x  =  b.  A  system  of  equations  then  result  of  the  form  (a  two-step 
process) : 


*  n 

(Step  1)  A  *  T  =  B  *  T 

x,y  x,y 

n+1  * 

(Step  2)  C  *  T  =  D  *  T 

x,y  x,y 


where  A,  B,  C,  and  D  are  all  tridiagonal  matrices. 

*/ 

#include  "fluid. h" 
double  adiO 

int  i,  j,  m,  n; 
double  max,  diff; 

/*  Next,  bring  in  the  values  of  u  (psi_y)  and  v  (-  psi_x)  from  the 
potential  flow  due  to  a  wall  jet.  Input  the  velocities  from  file 
main.c.  */ 

/*  STEP  1:  Solve  the  tridiagonal  system 
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*  *  *  n  n  n 

cl  T  +  c2  T  +  c3  T  =  c4  T  +  c5  T  +  c6  T 

x+l,y  x,y  x-l,y  x,y+l  x,y  x,y-l 

where  the  coefficients  cl  through  c6  are  determined  by  */ 

for  (i  =  0;  i  <=  nx  ;  i++){ 
for  (j  =  0;  j  <=  ny  ;  j++){ 

cl[i][j]  =  (xcon2  *  u[i][j])  -  (xconl  /  2.0); 
c3[i][j]  =  (-xcon2  *  u[i][j])  -  (xconl  /  2.0); 
c4[i]  [j]  =  (-ycon2  *  v[i][j])  +  (yconl  /  2.0); 
c6[i][j]  =  (ycon2  *  v[i][j])  +  (yconl  /  2.0); 

/*  printf ("*/,d  7.d  7,1  Af  %1.2f  7.1. 4f  %1.4f  7.1. 2f  7.1.4f\n",  i,  j, 
cl  [i]  [j]  ,  c2,  c3  [i]  [j]  ,  c4[i][j],  c5,  c6[i][j]);  */ 

> 

> 

/*  c2  and  c5  are  defined  in  initialization  file.  */ 

/*  Start  with  the  corner.  */ 

j  =  0; 
i  ■  1; 

tempaCi-1]  =  yconl  *  temp[i]  [j+1]  +  c5  *  temp[i]  [j]  -  (yconl  * 
hy  *  f[i])  -  c3[i][j]  *  temp[i-l]  [j] ; 

A[l]  [i-1]  =  c2; 

A [2]  [i-1]  =  cl  [i]  [j]  ; 

/*  printf ("7.d  7.d  7.1.4f\n",  i,  j,  tempa[i-l]);  */ 

/*  Next,  go  along  the  surface  to  the  other  corner.  */ 
for  (i=2;i<nx-l;  i++){ 

tempa[i-l]  =  yconl  *  temp [i]  [j+1]  +  c5  *  temp[i]  [j]  -  (yconl  * 
hy  *  f  [i] ) ; 

A[0]  [i-1]  =  c3[i]  [j]  ; 

A[l]  [i-1]  =  c2; 

A  [2]  [i-1]  =  cl  [i]  [j]  ; 

/*  printf ("7.d  7.d  7»1.4f\n",  i,  j,  tempa[i-l]);  */ 

> 

i  =  nx  -  1; 

tempa[i-l]  =  yconl  *  temp[i][j+l]  +  c5  *  temp[i]  [j]  -  (yconl  * 
hy  *  f[i])  -  cl  [i]  [j]  *  temp[i+l]  [j]  ; 
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A [0]  [i-1]  =  c3[i]  [j]  ; 

A [1] [i-1]  =  c2; 

/*  printf("70d  7.d  %1.4f\n",  i,  j,  tempa[i-l]);  */ 

/*  for  (m  =  0;  m  <  3;  m++H 

for  (n  =  0;  n  <  ny  -  1;  n++){ 

printf("7,d  7.d  7,1.4f\n",  m,  n,  A[m][n]); 

> 

}  */ 

size  =  maxxvec; 

/*  Use  tridiag.c  to  solve  A  T*  =  Tn,  where  A  =  is  tridiagonal, 

Tn  is  the  known  temperature  (tempa) ,  and  T*  is  the  step  temperature .  */ 

tridiag(A,  tempa,  size); 

for  (i  =  1;  i  <=  size  ;  i++){ 
tstar[i][j]  =  xx [i-1]  ; 

/*  printf  ("7.d  7.d  7.1.4f\n",  i,  j,  tstar[i]  [j] ) ;  */ 

} 

/*  This  is  the  solution  to  the  surface  row.  The  tstar  matrix  is 
updated  with  the  new  surface  temperature.  Now,  work  on  the  inside 
of  the  regime.  Continue  to  update  the  tstar  matrix  as  each  row  is 
solved.  */ 

for  (j  =  1;  j  <  ny;  j++M 
i  =  l; 

tempb [i-1]  =  c4[i]  [j]  *  temp[i]  [j+1]  +  c5  *  temp[i][j]  +  c6[i]  [j] 

*  temp[i]  [j-1]  -  c3[i]  [j]  *  temp[i-l]  [j]  ; 

A[l]  [i-1]  =  c2; 

A  [2]  [i-1]  =  cl  [i]  [j]  ; 

/*  printf  ("7.d  %d  %1.4f\n" ,  i,  j,  tempb  [i-1]);  */ 
for  (i=2;i<nx-l;  i++){ 

tempb  [i-1]  =  c4[i]  [j]  *  temp  [i]  [j+1]  +  c5  *  temp[i]  [j]  +  c6[i]  [j] 

*  temp [i]  [j-1]  ; 

A[0]  [i-1]  =  c3[i]  [j]  ; 

A[l]  [i-1]  =  c2; 

A  [2]  [i-1]  =  cl  [i]  [j]  ; 

> 
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i  =  nx  -  1; 

tempb[i-l]  =  c4[i]  [j]  *  temp[i]  [j+1]  +  c5  *  temp[i][j]  +  c6[i]  [j] 

*  temp[i][j-l]  -  cl[i][j]  *  temp[i+l]  [j]  ; 

A  [0]  [i-1]  =  c3[i]  [j]  ; 

A  Cl] [i~l]  =  c2; 

/*  printf ("7,d  %d  7.1.4f\n",  i,  j,  tempb[i-l]);  */ 

/*  Use  crout(A)  and  croutslvCL,  U,  tempb)  to  solve  for  row  j.  Then 
update  the  temp  matrix.  */ 

tridiag(A,  tempb,  size); 

for  (i  =  1;  i  <=  size  ;  i++){ 
tstar[i][j]  =  xx[i-l]; 

/*  printf  ("%d  7.d  7,1.4f\n",  i,  j,  tstar  [i]  [j]  )  ;  */ 

y 

> 

/*  Finally,  solve  along  the  bottom  row.  */ 

j  =  ny; 
i  =  1; 

tempo  [i-1]  =  c5  *  temp[i]  [j]  +  yconl  *  temp  [i]  [j-1]  -  c3[i]  [j] 

*  temp  [i-1]  [j]  ; 

A[l]  [i-1]  =  c2; 

A  [2]  [i-1]  =  cl[i]  [j]  ; 

/*  printf  ("*/,d  */,d  */,1.4f\n",  i,  j,  tempo  [i-1]);  */ 

for  (i  =  2;  i  <  nx  -  1;  i++){ 

tempo  [i-1]  =  c5  *  temp[i]  [j]  +  yconl  *  temp [i]  [j-1]  ; 

A[0]  [i-1]  =  c3[i]  [j]  ; 

A[l]  [i-1]  =  c2; 

A  [2]  [i-1]  =  cl  [i]  [j]  ; 

/*  printf  ("7,d  °/,d  7,1.4f\n",  i,  j,  tempo  [i-1]);  */ 

y 

i  =  nx  -  1; 

tempc[i-l]  =  c5  *  temp[i]  [j]  +  yconl  *  temp [i]  [j-1]  -  cl[i][j] 

*  temp  [i+1]  [j]  ; 

A[0]  [i-1]  =  c3[i]  [j]  ; 

A[l]  [i-1]  =  c2; 

/*  printf  ( "7, d  7.d  7.1.4f\n",  i,  j,  tempo  [i-1]);  */ 
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for  (m  =  0;  m  <  3;  m++){ 

for  (n  =  0;  n  <  ny  -  1;  n++){ 

/*  printf  ("7.d  %d  %1.4f\n",  m,  n,  A [m]  [n]  )  ;  */ 

> 

> 

for  (i  =  0;  i  <  nx  -  1;  i++){ 

/*  printf  ('7.1.4f\n\  tempo [i]);  */ 

> 

tridiag(A,  tempo,  size); 

for  (i  =  1;  i  <=  size  ;  i++){ 
tstar[i]  [j]  =  xx[i-l]  ; 

/*  printf  ("7.d  7.d  %1.4f\n",  i,  j,  tstar  [i]  [j] )  ;  */ 

y 

/*  Print  the  results  of  the  first  step  of  tstar.  */ 

/*  printf  C"x[i]  y[j]  tstar [i]  [j]\n")  ;  */ 
for  (i  =  0;  i  <=  nx;  i++){ 
for  (j  ■  0;  j  <=  ny;  j++){ 

/*  printf  ("*/,1.2f  7.1. 2f  7.1.4f\n",  x[i],  y[j]  ,  tstar  [i]  [j] )  ;  */ 

} 

} 

/*  STEP  2:  Repeat  in  the  other  direction.  Now  solve 

n+1  n+1  n+1  *  *  * 

dl  T  +  d2  T  +  d3  T  =  d4  T  +  d5  T  +  d6  T 

x,y+l  x,y  x,y-l  x+l,y  x,y  x-l,y 

where  the  coefficients  dl  through  d6  are  determined  by  */ 

for  (i  =  0;  i  <=  nx  ;  i++){ 
for  (j  =  0;  j<=  ny  ;  j++){ 

dl [i]  [j]  =  (ycon2  *  v[i][j])  -  (yconl  /  2.0); 
d3[i][j]  =  (-ycon2  *  v[i][j])  -  (yconl  /  2.0); 
d4[i][j]  =  (-xcon2  *  u[i]  [j])  +  (xconl  /  2.0); 
d6[i][j]  =  (xcon2  *  u[i][j])  +  (xconl  /  2.0); 

/*  printf ("%d  7.d  7.1. 4f  7.1- 4f  7.1. 4f  7.1- 4f  7.1. 4f  7.1.4f\n\  i,  j, 
dl  [i]  [j]  ,  d2,  d3[i]  [j]  ,  d4[i][j],  d5,  d6[i][j]);  */ 

} 
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> 

/*  d2  and  d5  are  defined  in  initialization  file.  */ 

/*  Only  one  loop  is  needed,  as  the  equations  axe  the  same  for  all 
i  columns.  */ 

/*  printf  C"x[i]  y[j]  temp[i]  [j]\n")  ;  */ 
size  =  maxx; 
max  =  0. ; 

for  (i  =  1;  i  <  nx;  i++){ 

j  =  0; 

tstarn[j]  =  d4[i][j]  *  tstar[i+l]  [j]  +  d5  *  tstar[i]  [j]  +  d6[i][j] 

*  tstar[i-l]  [j]  -  (yconl  *  hy  *  f[i]); 

C[l]  [0]  =  d2; 

C  [2]  [0]  =  -yconl ; 

for  (j  =  1;  j  <  ny;  j++){ 

tstarn[j]  =  d4[i]  [j]  *  tstar[i+l]  [j]  +  d5  *  tstar[i]  [j]  +  d6[i]  [j] 

*  tstar  [i-1]  [j]  ; 

CEO]  Cj]  =  d3[i]  [j]  ; 

C[l]  [j]  =  d2; 

C [2]  [j]  =  dl  [i]  [j]  ; 

> 

j  =  ny; 

tstarn[j]  =  d4[i]  [j]  *  tstar [i+1]  [j]  +  d5  *  tstar [i]  [j]  +  d6[i][j] 

*  tstar  [i-1]  [j]  ; 

C [0] [ny]  =  -yconl ; 

C[l]  [ny]  =  d2; 

tridiag(C,  tstarn,  size); 

for  (j  =  0;  j  <  size  ;  j++){ 

diff  =  fabsC  temp[i][j]  -  xx[j]); 
if  (diff  >  max)  { 
max  =  diff; 

/*  printf  ("new  max  [*/,d,7,d]  =  7,f\n"  ,i,  j  ,max)  ;  */ 

> 

temp  [i]  [j]  =  xx [j]  ; 

> 

> 

/*  Print  the  results  for  the  second  step  and  begin  the  next  iteration 


103 


with  the  new  temp.  */ 


for  (i  =  0;  i  <=  nx;  i++){ 
for  (j  =  0;  j  <=  ny;  j++){ 

/*  printf ("*/,1.2f  */,1.2f  */,1.4f\n",  x[i],  y[j],  temp[i][j]>;  */ 

> 

> 


return (max) ; 


/*  this  file:  tridiag.c  Solve  tridiagonal  system  A  x  =  b.  */ 

#include  "fluid.h" 

tridiag(t,  b,  size) 

double  t [3] [maxx] ,  b [maxx] ; 
int  size; 

{ 

double  r; 

for  (j  =  1;  j<  size;  j++){ 
r  =  t[0][j]  /  t  [1]  [j-l] ; 
t  Cl]  [j]  -=  r  *  t  [2]  [j-l]  ; 
b[j]  -=  r  *  b[j-l]  ; 

> 

b  [size-1]  /=  t[l]  [size-1]  ; 

/*  printf("7,d  %1.4f\n",  size-1,  b[size-l]);  */ 

for  (j  =  size-2;  j  >=  0;  j — ){ 

b[j]  =  (b[j]  -  (t [2]  [j]  *  b[j+l]))  /  t [1]  [j]  ; 

/*  printf ("7,d  7,1.4f\n",  j,  b[j]);  */ 

> 

/*  Solution  is  stored  in  xx[] .  */ 

for  (j  =  0;  j<  size;  j++){ 
xx [j]  =  b[j]  ; 

> 

> 
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/*  this  file:  surf ace. c  This  is  a  file  to  determine  the  stream 

function  from  the  surface  boundary  layer.  */ 


#include  "fluid. h" 
surface () 

for  (i  =  0;  i  <  nx;  i++){ 

delta[i]  =  pow(eps ,-l . 0/3.0)  *  pow (-temp [i] [0]  , 1 .0/3.0)  * 
pow (Tx [i] , 1 . 0/3 . 0) ; 

seta[i]  =  -  pow (eps, -2 .0/3.0)  *  pow (-temp [i] [0] ,2. 0/3.0)  * 
pow(Tx[i] ,-l. 0/3.0) ; 

for  (j  =  0;  j  <=  ny;  j++){ 

spsi[i][j]  =  seta[i]  *  (  (1. 0/2.0)  *  (1.0  - 
exp(-  pow(y [j] /(Pr*delta[i] ) ,2.0)))  +  ((sqrt (3. 14159265) /2.0) 

*  (y  [j]/(Pr*delta[i] ))  *  erfc(y[j]/(Pr*delta[i])))  ); 

/*  printf ("%1.3f  %1.3f  %1.6f\n",  x[i],  y[j],  spsi[i][j]);  */ 

y 

y 

y 
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/*  this  file:  wall.c  This  is  a  file  to  determine  the  stream 
function  in  the  wall  boundary  layer.  */ 

#include  "fluid. h" 

wall() 

double  coni,  con2,  psicon,  Delta,  feta; 
double  getgO; 

coni  =  (5. 0/2.0)  *  F; 

con2  =  2.0  /  pow( coni, 1. 0/4.0) ; 

for  (j  =  1;  j  <=  ny;  j++){ 

Delta  =  con2  *  pow(y [j] ,3. 0/4.0) ; 
psicon  =  -  2.0  *  pow((conl*y[j]) ,1.0/4. 0) ; 

for  (i  =  0;  i  <=  nx;  i++){ 
eta  =  x[i]  /  (Delta  *  Pr) ; 

/*  printf ("eta  =  */.1.3f\n",  eta);  */ 
g  =  getg(eta) ; 
feta  =  g  *  g; 

/*  printf  ("feta  =  */.1.6f\n",  feta);  */ 
wpsi[i][j]  =  psicon  *  feta; 

/*  printf  ( "7.1. 3f  %1.3f  7.1.6f\n",  x[i]  ,  y[j]  ,  wpsi[i][j]);  */ 

> 

> 

> 
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/*  This  file:  getg.c  -  This  is  a  file  to  determine  eta  using 
Newton's  Method,  to  be  used  in  wall.c  */ 

#include  "fluid. h" 

double  getg(eta) 

double  eta; 

double  gO,  gprime,  etaO,  tol; 

tol  =  0.000001; 

if  (eta  >  7.8)  return(l.O); 

if  (eta  <  exp (1.0)) 
gO  =  eta  /  3.0; 
else 

{gO  =  1.0  -  sqrt(3.0  /  exp(2.0  *  eta  -  3.14159265  /  sqrt (3.0))) ;} 
/*  printfC'gO  =  %1.6f\n",  gO) ;  */ 


g  =  gO; 
do{ 

gO  =  g; 

gprime  =  (1  -  g  *  g  *  g)/  3.0; 

etaO  =  Iog((sqrt(1.0+g0+g0*g0))/(1.0-g0))  +  (sqrt(3.0))  * 
atan(((sqrt(3.0))/3.0)*(l .0+2.0*g0))  -  (sqrt (3.0))* (3. 14159265/6.0) ; 

g  =  gO  +  (eta  -  etaO)  *  gprime; 

/*  printf  ("g  =  */.1.6f\n",  g) ;  */ 

}  while  (fabs(g  -  gO)  >  tol); 
return (g) ; 


> 
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/*  this  file:  main.c  -  main  program  for  uniform  stream  function.  */ 


#def ine  MAIN 
#include  "fluid. h" 

#def ine  LINESIZE  128 

#def ine  eps  (sqrt(3.14159265)*(1.0-sqrt(2.0)/2.0)/2.0) 
main() 

-c 

psiinit () ; 
tx(); 

strength () ; 

jet() ; 

green C); 
wallO; 
surf ace () ; 

for  (i  =  0;  i  <=  nx;  i++){ 
for  (j  =  0;  j  <=  ny;  j++){ 

unifpsi  [i]  [j]  =  jpsi[i][j]  +  gpsi[i][j]  +  wpsi[i][j] 

+  spsi  [i]  [j]  -  jpsi  [0]  [j]  -  gpsi  [i]  [0]  ; 

printf  ("y,1.3f  */.1.3f  */,1.6f  %1.6f  */.1.6f  */,1.6f  */,1.6f\n", 
x[i],  y[j] ,  jpsi  [i]  [j] ,  gpsi[i][j],  spsi[i][j],  wpsi[i][j], 
unifpsi  [i]  [j]); 

> 

> 

> 
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/*  this  file:  psiinit.c  -  initialization  for  uniform  stream  function.  */ 


#include  "fluid. h" 
#def ine  LINESIZE  128 

psiinitO  { 

FILE  *f ile; 

char  line [LINESIZE] ; 

nx  =  ny  =  maxx  -  1; 
hx  =  hy  =  3.0  /  nx; 


for  (i  =  0;  i  <=  nx  ;  i++){ 
x[i]  =  y [13  =  i  *  hx; 

> 

file  =  f open ("input.txt",  "r"); 

for  (i  =  0;  i  <=  nx;  i++){ 
forCj  =  0;  j  <=  ny;  j++){ 

fscanf  (f  ile,  "°/.lf*/.lf7,lf\n"  ,  &Cx[i]),  &(y[j]),  &(temp[i]  [j]  ))  ; 

> 

> 

/*  printf  C'"/,1.3f  7,1. 3f  %1.6f\n",  x[l],  y[2],  temp  [1]  [2]);  */ 

} 
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/*  this  f ile:unitemp. c  -  Calculates  the  uniform  temperature  solution 
by  matching  Tinner  +  Touter  -  Tmatch.  */ 

#define  MAIN 
#include  "fluid.h" 

#def ine  LINESIZE  128 

#define  eps  (sqrt (3. 14159265) *(1.0-sqrt (2. 0)/2.0)/2.0) 

mainC) 

{ 

double  unitemp  [maxx]  [maxx]  ,  teta; 

nx  =  ny  =  maxx  -  1; 
hx  =  5.0  /  nx; 

for  (i  =  0;  i  <=  nx  ;  i++){ 
x [i]  =  y [i]  =  i  *  hx; 

> 

for  (i  =  0;  i  <=  nx  ;  i++){ 
for (j  =  0;  j  <=  ny;  j++){ 

scanf  ("7.*lf  7.*lf  7.1f  "  ,  &(temp  [i]  [j]  )  )  ; 

> 

} 

for  (i  =  0;  i  <=  nx;  i++){ 
if  (i  ==  0) 

Tx[i]  =  (4.0*temp[l]  [0]  -3 . 0*temp [0]  [0] -temp [2]  [0]  )/(2 .0*x  [i+1]  )  ; 
else  if  (i  ==  nx) 

Tx[i]  =  (temp  [nx]  [0] -temp  [nx-1]  [0])  /  (x[i]-x[i-l]  )  ; 
else 

Tx[i]  =  (temp [i+1]  [0] -temp [i— 1]  [0])  /  (2.0*  (x[i+l] -x [i]  ))  ; 

/*  printf ("7.1 .2f  7.1. 4f  7.1.4f\n",  x[i],  temp[i][0],  Tx[i]);  */ 

> 

for  (i  =  0;  i  <=  nx  ;  i++){ 
f °r(j  =  0;  j  <=  ny;  j++){ 

teta  =  y [j]  /  (Pr  *  pow ( (-temp [i]  [0] )/ (eps  *  Tx[i]  *  Tx[i]), 

1. 0/3.0)  ); 

/*  printf ("7.1. 2f  7.1. 2f  7,1.15f\n",  x[i],  y[j],  teta);  */ 

unitemp [i]  [j]  =  temp[i]  [j]  +  (Pr  /  (4.0*eps))  *  temp[i][0]  * 

(  ((sqrt (3. 14159265))  *  (erf c (teta)))  *  (0.5  +  teta  *  teta)  - 
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(teta  *  exp(-  teta  *  teta))  ); 

printf  ("%1.2f  %1.2f  %1.5f\n",  x[i],  y[j],  unitemp  [i]  [j]  )  ; 

> 

> 

> 
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